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Abstract 

The  numerical  approximation  of  definite  integrals,  or  quadrature,  often  involves 
the  construction  of  an  interpolant  of  the  integrand  and  subsequent  integration  of  the 
interpolant.  It  is  natural  to  rely  on  polynomial  interpolants  in  the  case  of  one  di¬ 
mension;  however,  extension  of  integration  of  polynomial  interpolants  to  two  or  more 
dimensions  can  be  costly  and  unstable.  A  method  for  computing  surface  integrals 
on  the  sphere  is  detailed  in  the  literature  (Reeger  and  Fornberg,  Studies  in  Applied 
Mathematics ,  2016).  The  method  uses  local  radial  basis  function  (RBF)  interpola¬ 
tion  to  reduce  computational  complexity  when  generating  quadrature  weights  for  a 
particular  node  set.  This  thesis  expands  upon  the  same  spherical  quadrature  method 
and  applies  it  to  an  arbitrary  smooth  closed  surface  defined  by  a  set  of  quadrature 
nodes  and  triangulation. 

Key  words:  Quadrature,  Radial  Basis  Functions,  RBF,  Numerical  Integration, 
Smooth  Surface,  Interpolation,  Closed  Surface 
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Radial  Basis  Function  Based  Quadrature  Over 
Smooth  Surfaces 

I.  Literature  Review 

1.1  Background  Motivation 

An  increasing  number  of  applications,  arising  for  example  in  geophysics,  require 
partial  differential  equations  (PDEs)  to  be  solved  on  irregular  geometries  [1].  A 
numerical  solution  of  PDEs  is  often  paired  with  numerical  quadrature  over  smooth 
(closed)  surfaces  in  order  to  obtain  integrated  quantities,  such  as  total  energy,  average 
temperature,  etc.  [2],  Accurate  approximation  of  PDEs  and  integrals  on  a  surface 
often  requires  the  construction  of  node  sets  featuring  spatially  varying  densities  in 
order  to  capture  rapidly  changing  features  (viz.  regions  of  large  curvature).  This 
thesis  presents  an  approach  for  calculating  quadrature  weights  for  such  node  sets  to 
be  used  in  integral  approximations  over  a  smooth  closed  surface  S.  That  is, 

Mf)  -  ///(*.» (1) 

S 

where  f(x,y,z)  is  a  scalar  function  known  at  some  N  quadrature  nodes  and  WV: 
are  the  corresponding  quadrature  weights. 

Earlier  literature  focused  on  computing  surface  integrals  on  the  sphere,  or  spherical 
quadrature,  over  very  specific  node  sets,  partnered  with  tabulated  weights  for  select 
values  of  N  (the  total  number  of  nodes)  [3,  4,  5].  More  recent  literature  borrowed 
concepts  from  radial-basis-function-generated  finite  differences  (RBF-FD)  applied  to 
spatially  variable  node  sets  for  spherical  quadrature  [2],  This  more  recent  technique 
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for  spherical  quadrature  is  expanded  upon  in  Chapter  II  and  adapted  for  use  in 
approximating  integrals  of  arbitrary  smooth  closed  surfaces.  The  need  for  approxi¬ 
mations  to  surface  integrals  on  surfaces  other  than  the  sphere  arises  in  applications 
that  vary  from  geophysics,  where  the  earth  may  be  treated  as  an  ellipsoid  rather  than 
a  sphere,  to  biology,  where  particle  interactions  may  require  models  more  complicated 
than  colliding  spheres  (see,  for  example,  [6]  and  references  therein). 

This  chapter  introduces  the  inherent  benefits  of  using  radial  basis  functions  (RBFs) 
as  an  approximation  basis  for  node  sets  featuring  spatially  variable  density.  These 
benefits,  which  include  high  orders  of  accuracy  when  integrating  the  interpolant  (e.g. 
0(N~3'5)  for  a  sphere  on  nearly  uniformly  spaced  points)  and  relatively  cheap  com¬ 
putation  cost  (e.g.  O(NlogN)  for  the  techniques  borrowed  from  [2]),  are  why  RBFs 
were  chosen  for  use  during  interpolation  in  this  thesis. 

1.2  Radial  Basis  Functions 

The  approximation  of  a  multivariate  function  /(x),  x  e  (d  the  dimension), 
can  be  accomplished  in  many  ways  -  for  example,  the  use  of  multivariate  polynomials 
as  a  basis  [7].  When  moving  beyond  one  dimension,  Mairhuber  [8]  showed  that 
interpolation  with  a  set  of  basis  functions  that  do  not  depend  on  the  data  locations  is 
an  ill-posed  problem.  That  is,  an  infinite  number  of  configurations  of  the  points  where 
the  interpolation  conditions  are  enforced  lead  to  a  singular  linear  system  of  equations 
on  the  interpolation  coefficients.  Multivariate  polynomial  interpolation  suffers  from 
this  problem.  Radial  basis  function  interpolation,  on  the  other  hand,  uses  (as  the  set 
of  basis  functions)  translates  of  a  single  function  {4>(ri)}f=1,  where  r y  =  ||x  —  Xj||  is 
the  distance  of  any  point  x  in  the  observational  domain  to  the  point  x,:  (where  an 
interpolation  condition  will  be  enforced).  As  such,  RBFs  do  not  suffer  from  the  same 
problem  as  multivariate  polynomials  and  become  useful  in  this  regard. 
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Not  only  can  radial  basis  functions  be  applied  for  any  number  of  dimensions,  but 
they  also  have  excellent  approximation  properties.  They  can  be  applied  on  meshless 
data  sets,  allowing  a  user  to  easily  enhance  interpolation  properties  over  certain  data 
features.  For  example,  there  may  be  unique  results  occurring  for  a  certain  region 
of  interest.  Clustering  data  nodes  near  the  region  of  interest  reduces  uniformity  in 
data  fidelity.  However,  radial  basis  functions  have  the  ability  to  work  regardless  of 
this  constraint.  Examples  of  RBFs  can  be  found  in  table  1  [9,  10],  where  £  is  a 
non-negative  parameter  for  the  infinitely  smooth  RBFs  that  affects  the  RBF  shape. 

Table  1.  Types  of  Radial  Basis  Functions  </>(r) 


Piecewise  Smooth  (Conditionally  Positive  Definite) 


MN 

Monomial 

|  j,  |2m+l 

TPS 

thin  plate  spline 

r  |2mln  r 

Infinitely  Smooth  (Positive  Definite) 

MQ 

Multiquadric 

a/1  +  (er)2 

IQ 

Inverse  quadric 

i 

1 +(er)2 

IMQ 

Inverse  MQ 

1 

y/l+(er)2 

GA 

Gaussian 

e-(-)2 

BE 

Bessel 

Jd  Aer) 

d  -i 

CP2 


1.3  History  of  RBFs 

Radial  basis  functions  are  a  more  recent  development  in  mathematics.  Uncon¬ 
ditional  nonsingularity  of  the  interpolation  problem  for  a  specific  case  was  given  by 
German  scientist  Bochner  in  1933  [11].  This  was  followed  by  the  use  of  positive  def¬ 
inite  functions,  a  class  of  functions  to  which  many  RBFs  belong,  in  specific  metric 
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spaces  in  1938  by  Schoenberg  [12].  It  was  not  until  1971,  however,  that  Hardy  applied 
multiquadric  (MQ)  functions  to  applications  in  topography  given  coordinate  data  for 
irregular  surfaces  [13].  The  multiquadric  functions  are  radial  functions  of  the  form 
indicated  in  table  1. 

In  Hardy’s  work,  RBFs  were  used  to  interpolate  irregular  smooth  surfaces  given 
discrete  data  [13].  In  1990,  he  published  another  article  summarizing  MQ  and 
multiquadric-biharmonic  (MQ-B)  methods  for  use  in  such  applications  as  geodesy, 
geophysics,  surveying  and  mapping,  geography,  remote  sensing,  signal  processing, 
and  more  [1].  After  this  emergence  of  MQ  equations,  radial  basis  functions  became 
increasingly  popular. 

Also  in  1990,  Kansa  published  a  paper  using  multiquadrics  as  a  basis  for  approx¬ 
imating  the  solution  to  partial  differential  equations  [14].  He  pointed  out  that  these 
functions  are  continuously  differentiable  and  integrable  over  a  meshless  grid,  and  are 
therefore  useful  for  representing  other  functions  in  the  setting  of  integral  and  differ¬ 
ential  equations.  Their  simple  extension  to  higher  dimensions  makes  multiquadrics 
and  other  radial  basis  functions  great  candidates  for  solving  PDEs. 

The  construction  of  a  RBF  interpolant  can  be  computationally  intensive,  par¬ 
ticularly  if  the  node  set  S n:  {xj}[^1  (where  the  interpolation  conditions  are  to  be 
enforced)  is  large.  In  the  presence  of  N  interpolation  conditions  on  N  nodes,  a  dense 
linear  system  of  equations  of  size  IV  x  N  must  be  solved.  Further  complicating  the 
matter  is  that  these  interpolants  are  often  used  to  construct  approximate  differential 
operators  (as  in  the  case  of  [14])  whose  application  is  through  a  dense  matrix  multipli¬ 
cation.  This  computational  complexity  inspired  the  development  of  RBF-generated 
finite  differences  (RBF-FD)  [15,  16,  17,  18],  which  considers  for  each  node  Xj  €  S n  a 
subset  Mn  C  of  n  points  nearest  x,:  and  constructs  an  interpolant  (and  subsequent 
approximate  differential  operator)  on  that  set.  This  reduces  computational  cost  to  N 
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linear  systems  of  size  n  x  n  and  an  approximate  differential  operator  that  is  sparse 
[19].  The  differences  between  global  RBF  approximations  and  RBF-FD  has  analogy 
to  those  between  pseudo-spectral  methods  and  finite  differences  for  approximating 
derivatives. 

1.4  Numerical  Quadrature  over  Surfaces 

The  strong  approximation  properties  of  radial  basis  functions  in  interpolation 
and  the  solution  of  PDEs  led  to  their  consideration  for  approximations  to  surface 
integrals.  Reeger  and  Fornbeg  discuss  numerical  quadrature  over  the  surface  of  a 
sphere  [2],  forming  the  building  blocks  for  this  thesis.  As  a  result,  much  of  their  work  is 
reiterated  in  Chapter  II  for  completeness.  Various  integration  methods  on  the  sphere 
are  written  by  [5,  20,  21,  22,  23].  For  instance,  [23]  proposes  spherical  harmonics  as 
a  set  of  basis  functions  for  interpolation  and  analyzes  the  variability  in  results  for 
using  different  quasi-uniform  node  sets  on  §2.  On  the  other  hand,  when  applying  the 
associated  quadrature  weights  to  approximate  a  surface  integral,  Womersley  and  Sloan 
[5]  suggest  the  use  of  exact  cubature  weights  (quadrature  weights  in  M3)  associated 
with  polynomials  of  certain  degrees,  while  also  analyzing  the  use  of  various  quasi¬ 
uniform  node  sets  on  §>2. 

Some  of  the  methods  for  numerically  computing  surface  integrals  on  the  sphere 
have  been  adapted  to  other  surfaces  [24,  25].  Just  as  this  thesis  adapts  from  a  spher¬ 
ical  quadrature  method  [2]  to  approximate  an  integral  over  a  smooth  closed  surface 
S  C  M3,  Atkinson  [24]  also  adapts  spherical  quadrature  methods  to  solve  integral 
equations  defined  over  simple  smooth  surfaces.  He  presents  two  general  methods: 
product  Gaussian  quadrature  and  finite  element  integration.  The  product  Gaussian 
quadrature  uses  Gauss-Legendre  nodes  and  quadrature  weights  for  the  unit  sphere 
§2,  with  convergence  rate  0((2A^1//2  —  1)_“)  on  nearly  uniform  node  sets  ( u  the  nuni- 
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ber  of  times  /  in  (1)  is  continuously  differentiable).  The  finite  element  integration 
discretizes  §2  into  triangles,  with  convergence  rate  0(1F_1)  (K  the  number  of  trian¬ 
gles)  when  using  spherical  triangles,  and  0(K~2)  when  mapping  spherical  triangles 
to  planar  triangles.  He  interpolates  integrals  over  simple  smooth  surfaces  using  a 
one-to-one  mapping  with  §2  coupled  with  one  of  his  two  methods. 

Fuselicr  et  al.  [25]  developed  quadrature  methods  that  extend  to  manifolds  M 
that  are  either  homogeneous  (including  §2)  or  diffeomorphic  to  homogeneous  spaces. 
Their  methods  result  in  a  convergence  rate  of  roughly  0(N~2)  for  sufficiently  smooth 
integrands. 

Additional  types  of  quadrature  for  surfaces  in  M3  can  be  found  at  [26,  27,  28]. 
For  example,  [26]  discusses  numerical  quadrature  for  piecewise  smooth  surfaces  using 
polynomial  interpolants,  while  [27]  couples  Thin-Plate  Spline  interpolation  (see  table 
1)  with  Green’s  integral  formula  [29]  for  a  meshless  cubature  in  M3.  Quadrature 
involving  a  Galerkin  discretization  of  a  boundary  integral  equation  with  a  weakly 
singular  kernel  is  adapted  from  electromagnetics  to  a  general  framework  in  [28] . 

In  any  case,  a  common  theme  throughout  the  development  of  each  of  these  meth¬ 
ods  is  that  they  are  fixed  to  particular  node  sets  featuring  near-uniform  density  in 
order  to  achieve  stability.  Some  of  these  can  achieve  high  orders  of  accuracy,  even 
spectral,  on  the  near-uniform  data  sets,  but  at  the  expense  of  a  lengthy,  often  in¬ 
tractable,  process  for  constructing  the  node  sets  themselves.  For  further  discussion 
of  node  sets  on  §>2,  see  section  2.1.1. 

1.5  Organization  of  Thesis 

The  aim  of  this  thesis  is  to  evolve  the  spherical  quadrature  method  described  in 
the  literature  [2]  and  apply  it  to  approximate  the  surface  integral  of  a  scalar  function 
over  an  arbitrary  smooth  closed  surface.  Chapter  II  walks  through  the  steps  involved 
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in  applying  the  spherical  quadrature  method,  from  discretizing  a  sphere  surface  to 
computing  the  quadrature  weights.  This  is  followed  in  Chapter  III  by  a  similar 
procedure  applied  to  an  arbitrary  smooth  surface.  Chapters  IV  and  V  include  some 
results  of  the  methods  described  in  Chapter  III  and  possible  future  considerations. 
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II.  Spherical  Quadrature 


Much  of  the  basis  behind  the  method  for  numerical  quadrature  over  a  smooth 
closed  surface  comes  from  the  techniques  used  for  spherical  quadrature  discussed  in 
[2] .  The  purpose  of  this  chapter  is  to  explain  the  construction  of  the  set  of  weights  WV 
in  (1)  for  S  the  surface  of  a  sphere  with  radius  p  centered  at  the  origin.  The  N  nodes 
make  np  the  set  S]y\  {(aq,  ytl  Zi)}f=l.  This  notation  will  be  consistent  throughout  this 
chapter.  The  process  to  construct  the  quadrature  weights  can  be  summarized  by 

1.  Discretize  the  sphere  surface  using  a  node  set  Sn  and  triangulation  T. 

2.  For  each  triangle  T)c  G  T.  k  =  project  a  local  region  Snk  containing  n 

surface  nodes  (from  Sn)  and  corresponding  triangles  to  a  plane  tangent  to  the 
sphere. 

3.  Interpolate  the  transformed  integrand  over  n  projected  points  J\fn  C  Sn  in  the 
tangent  plane  using  radial  basis  functions  as  a  basis. 

4.  Integrate  the  interpolant  over  a  triangle  in  the  tangent  plane  region  to  get 
planar  quadrature  weights 

5.  Convert  the  planar  quadrature  weights  to  spherical  quadrature  weights  to  be 
used  in  (1) 

This  process  will  set  the  foundation  for  understanding  the  surface  quadrature  pre¬ 
sented  in  Chapter  III. 


2.1  Discretizing  the  Surface  S 


2.1.1  Types  of  Node  Sets. 

The  sphere  surface  §2  can  be  easily  parameterized,  and  hence  many  discretization 
processes  exist  to  yield  a  desired  set  of  surface  nodes  Sn  C  §2.  The  approximation 
in  (1)  requires  such  a  set  of  quadrature  nodes  Sjy  of  N  nodes  in  M3.  Various  types 
of  node  sets  exist  for  §2,  including  the  following  examples  illustrated  in  figure  1. 
Define  quasi-uniform  to  mean  that  the  spatial  density  of  a  set  of  points  on  §2  is  near 
uniform  (such  node  sets  are  necessary  when  considering  the  methods  discussed  in 
section  1.4).  That  is,  each  point  on  §>2  accounts  for  approximately  the  same  surface 
area.  Some  quasi-uniform  node  sets  include  Fibonacci,  Icosahedral,  minimum  energy, 
and  maximum  determinant.  The  Fibonacci  node  set  [30]  creates  a  mesh  such  that 
each  node  on  S>2  represents  roughly  the  same  area  (see  section  2.3.2  on  fill  distance). 
Meanwhile,  the  Icosahedral  node  set  discussed  in  [25]  subdivides  the  triangular  facets 
of  the  icosahedron.  Quasi-minimum  energy  nodes  in  [25],  on  the  other  hand,  treat 
the  nodes  as  point  charges  and  nearly  minimize  the  Riesz  energy  over  possible  node 
configurations  [31,  32],  Another  minimum  energy  node  set  locally  minimizes  potential 
energy  [5].  The  last  depicted  quasi-uniform  node  set,  the  maximum  determinant  node 
set  [5] ,  locally  maximizes  the  determinant  of  the  interpolation  matrix  on  the  basis  of 
spherical  harmonics.  These  node  distributions  vary  such  that  they  should  be  chosen 
based  on  the  application. 

To  illustrate  the  geometric  flexibility  of  the  present  spherical  quadrature  method, 
this  thesis  considers  the  Halton  node  set  that  maps  the  first  N  points  from  the  Halton 
sequence  in  M2  to  the  unit  sphere  §2  as  discussed  in  [23].  Likewise,  the  random  node 
set  maps  N  pseudo-random  (computer  generated)  points  from  the  standard  normal 
distribution  in  M2  to  §2  also  as  discussed  in  [23]. 

Each  of  these  node  sets  generates  a  different  set  of  quadrature  weights,  with  vary- 
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2.1.2  Delaunay  Triangulation. 

The  set  of  quadrature  nodes  S n  can  be  constructed  in  many  ways  as  seen  in  section 
2.1.1.  An  example  of  such  a  set  appears  in  figure  2a.  From  the  set  Sn  a  spherical 
triangulation  T  =  {rfc}f=1  is  constructed  [33]  such  that 

•  the  triangle  vertices  are  the  elements  of  Sn, 

•  the  triangle  edges  are  geodesics  between  the  vertices  of  the  triangles  (uniquely 
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defined  as  long  as  the  two  vertices  on  the  edge  are  not  separated,  in  angle,  by 
exactly  7 r), 

•  no  triangle  contains  an  element  of  Sjy  other  than  its  vertices, 

•  the  interiors  of  the  triangles  are  pairwise  disjoint, 

•  the  union  of  the  set  T  covers  S>2,  and 

•  no  circnmcircle  of  a  triangle  r*,  contains  an  element  of  S n  on  its  interior. 

An  example  Delaunay  triangulation  is  illustrated  in  figure  2b.  ft  is  a  generalization 
of  the  straight  line  dual  of  the  Voronoi  diagram,  formally  introduced  by  Delaunay  in 
1934  [34] .  The  particular  triangulation  listed  here  covers  the  convex  hull  of  the  sphere 
nodes,  and  is  further  discussed  by  Renka  in  [33].  Renka  also  includes  an  O(NlogN) 
algorithm  for  its  construction.  The  triangulation  can  also  be  generalized  to  higher 
dimensions  according  to  [35]. 

The  triangulation  T  allows  Is(f)  hr  (1)  to  be  written  as 

f(x,y,  z)dS 

Q  K=L  ^ 

such  that  each  r/c  €  T,  for  k  =  1 K,  can  be  considered  independently.  That  is, 
quadrature  can  be  performed  on  each  individual  triangle  before  adding  up  all  the 
quadrature  weights  to  yield  the  desired  weights  Wn  for  the  surface  integral  Xg{f). 

2.2  Projection  into  a  2-D  plane 

With  Sn  and  T  in  hand,  the  following  subsections  will  transform  a  three-dimensional 
region,  S^k,  on  a  sphere  surface  to  a  two-dimensional  region,  f4,  on  a  plane  tangent 
to  the  sphere. 


Zs (/)  =  II  f(x,y,z)dS  = 
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(a)  Nodes  Sjy  distributed  over  §2. 


(b)  Delaunay  triangulation  T  from 
section  2.1.2. 


(c)  Nearest  nodes  J\Tn  and  surface  region 
Snk  projected  to  tangent  plane. 


(d)  s(r /,  £)  interpolating  g(r /,£)  over 
planar  region  f 2^ 


Figure  2.  Progression  of  quadrature  steps  from  sections  2.1.2  through  2.2.3 
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2.2.1  Projecting  Neighboring  Nodes  to  a  Tangent  Plane. 

In  order  to  reduce  the  dimensionality  of  the  quadrature  problem  (1),  consider 
parameterizing  each  triangle  {r/,.}(L1  G  T  locally  through  a  change  of  variables.  The 
change  of  variables  is  performed  via  a  gnomonic  projection  [36]  into  a  plane  that  is 
tangent  to  the  sphere  at  a  point  in  each  triangle  Tk,  as  seen  in  figure  2c.  The  hallmark 
of  the  chosen  projection  must  be  that  geodesics  on  the  sphere  become  lines  in  the 
tangent  plane.  This  requirement  enforces  that  the  surface  area  of  each  triangle  in  T 
is  not  included  more  than  once,  and  that  the  surface  integral  over  a  spherical  triangle 
Tk  becomes  an  area  integral  over  a  planar  triangle  fk- 

To  form  the  projection  for  some  spherical  triangle  Tk ,  let  tabc  be  the  vertex 
representation  of  Tk  -  that  is,  x,i,xb,xc  G  Sn  are  the  vertices  of  Tk  -  and  define 

xm  =  [xMiVMi 

As  this  point  lies  below  the  sphere  surface,  the  spherical  triangle  midpoint  is 

~  r  ~  ~  ~  -i  ^-JVI 

Xm  :=  [XmjVm,  Zm\  =  y - 

IIxm||2 

(for  a  sphere  of  radius  p )  to  project  xm  to  the  sphere  surface  S.  The  equation  for  the 
tangent  plane  that  passes  through  xm  and  that  has  the  same  normal  vector  as  that 
of  the  sphere  surface  at  xm  is 

Xm(x  -  xm)  +  ym(y  -  ym)  +  zm(z  -  zm)  =  0 
since  Xm  —  0  is  the  surface  normal  at 

Consider  projecting  n  neighboring  nodes  (from  Sn  and  nearest  xm)  into  the  tan¬ 
gent  plane.  Call  this  set  of  neighboring  nodes  J\Tn.  Projecting  this  set  of  nodes  is  a 


zm]  — 


xa  +  xB  +  xc 
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necessity  when  approximating  a  surface  integral  over  an  entire  spherical  triangle 
via  an  area  integral  over  a  planar  triangle  ffc.  Let  Snk  represent  the  region  of  the 
sphere  surface  S  occupied  by  the  neighboring  nodes  in  A fn.  The  projection  occurs 
by  finding  the  intersection  between  the  plane  tangent  to  the  sphere  at  and  the 
line  passing  through  a  point  ( x,y,z )  G  S  near  (xm,Vm,Zm)-  That  is,  consider  the 
parameter  c  such  that 


xM(x  -  xM)  +  VM(y-  Vm)  +zm(z~  zm)  =  0 


-  - 

-  - 

X 

X 

y 

=  c 

y 

z 

z 

are  satished  simultaneously.  Therefore, 


x2m  +  Vm  +  =  P 2 _ 

x  Mx  +  yMy  +  zMz  xMx  +  yMy  +  zMz 


(2) 


Notice  that  a  singularity  in  c  occurs  where  Xm  ■  x  =  0,  for  •  the  vector  dot  product. 
That  is,  when  the  vector  in  the  direction  of  x  is  perpendicular  to  the  vector  in  the 
direction  of  the  triangle  midpoint  xM.  This  will  become  an  important  factor  when 
approximating  the  integral  of  /  in  the  tangent  plane  since  a  discrete  set  of  points 
from  Sn  near  r^,  but  not  within  r^,  will  need  to  be  projected  as  well.  As  a  result, 
only  points  Sn  within  a  90°  angle  of  5cm  are  considered  for  projection. 


2.2.2  Defining  a  Two-Dimensional  System. 

Now  that  the  region  Snk  about  a  spherical  triangle  Tk  G  T  has  been  projected  into 
a  tangent  plane,  it  is  natural  to  define  a  two-dimensional  coordinate  system  in  the 
plane  for  computation.  Let  represent  the  two-dimensional  planar  region  projected 
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and  transformed  from  Sak  ■  To  simplify  converting  the  three-dimensional  planar  coor¬ 
dinates  x:  (x,  y,  z )  to  two-dimensional  planar  coordinates  lj:  (r),  £),  consider  rotating 
the  coordinates  ( x,y,z )  so  that  Xm  is  located  along  the  z-axis  at  (0,  0,  p).  Further¬ 
more,  this  would  place  the  tangent  plane  at  parallel  to  the  (x,  y )  plane,  making 
all  5  set  equal  to  p.  This  can  be  done  through  a  series  of  rotation  matrices 


XM 

Vm 

0 

ZM 

n 

V*M+y2M 

V*M+y2M 

a/ 4+Vm 

p 

U 

p 

Vm 

xm 

0 

and  M2  — 

0 

i 

0 

a/  Xm+Vm 

a/  Xm+Vm 

0 

0 

1 

\J  Xm+Vm 

p 

0 

ZM 

P 

where  p  =  y/ x2M  +  fM  +  z2M .  So, 


M2Mi 


XMZM 

pV*M+y2M 

Vai 

y/*M+y2M 

XM 

P 


Vmzm 

P\f*M+yM 

XM 

VM 

P 


V^li+yh 

p 

0 


ZM 

P 


When  performing  the  rotation  on  the  triangle  midpoint, 


xm 

0 

m2m± 

Vm 

= 

0 

zm 

P 

verifies  the  intended  rotation.  Performing  the  transformation  from  three-dimensional 
coordinates  ( x,y,z )  G  Snk  to  two-dimensional  coordinates  (?/,£)  G  f4  (with  c  in  (2) 
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paired  with  the  rotation  matrices)  yields  the  transformation 


-  - 

-  - 

p[zM{xMX+yMy)-z{x2M+y2M)} 

V 

X 

y/  x2M+y2M(xMx+yMy+ZMz) 

£ 

=  cM2Mi 

y 

_ 

p2(xMy-yMx) 

\/  ^M+yM^MX+yMy+ZMz) 

7 

z 

P 

Just  as  intended,  the  last  element  of  the  planar  coordinates  is  the  radius  of  the 


sphere  and  can  therefore  be  dropped  from  the  two-dimensional  coordinate  system. 
This  yields  the  same  results  as  [2]: 


0 

1 


0 

0 


V 

£ 

7 


p[zM{xMx+yMy)-z(x2M+y2M)] 
yj  sm  +Vm  {xMX+yMy+zMz) 
_ p2(xMy-yMx) _ 

-  y/ sm  +Vm  {xMX+yMy+zMz)  . 


(3) 


for  x2M  +  yh  7^  0. 

That  is,  when  x«  is  already  at  (0,0,  p)  (or  (0,0,  — p))  and  before  applying  the 
rotation  matrices  M\  and  M2,  the  parameter  c  becomes 


_ f _ =f?_  =  P 

xMx  +  yMy  +  zMz  pz  z 
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This  results  in  (3)  becoming  (for  x2M  +  y2M  =  0)  [2] 


r  -| 

r 

1 

X 

V 

1 

0 

0 

<jl!  = 

0 

1 

0 

c 

y 

L  J 

L 

z 

px 

1 

0 

0 

z 

py 

0 

1 

0 

z 

_ 

_ 

pz 

z 

px 

z 

py 

z 


(4) 


Hence,  a  two-dimensional  coordinate  system  is  defined  for  each  tangent  plane  region 
Vtk  projected  from  Sak  about  rk  e  T  for  k  —  1, ...,  K. 


2.2.3  Radial  Basis  Function  Interpolation. 

With  a  two-dimensional  coordinate  system  defined,  it  is  imperative  to  understand 
radial  basis  function  (RBF)  interpolation  in  order  to  interpolate  the  integrand  g  in 

xSnk (/)  =  ff  f(x,y,z)dS  =  If  gfagjdrjdZ  =  Ink(g) 

for  Ink  the  area  integral  over  Qk  (see  figure  2d).  As  mentioned  in  section  1.3,  RBFs 
are  a  type  of  basis  function  that  can  be  used  in  interpolation  in  more  than  one  dimen¬ 
sion.  Given  a  set  of  points  (nodes)  A fn:  {xj}”=i  in  a  domain  and  corresponding 
measurements  or  samplings  of  an  unknown  scalar  function  g,  {g(xj)}”=1,  interpolation 
attempts  to  approximate  g  at  other  locations  (besides  those  from  A fn).  In  practice, 
A fn  will  consist  of  the  n  nodes  in  S n  nearest  xM  projected  to  the  plane  tangent  to  the 
sphere  at  Xm  and  transformed  to  the  two-dimensional  coordinate  system  (■ r] ,  £).  The 
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interpolation  is  done  by  forming  a  linear  combination  of  basis  functions  such  that  the 
interpolation  formula  (interpolant)  has  the  form 

n 

s(x)  =  ’ x  e  Rd  (5) 

j= 1 

Here,  d  is  the  dimension,  {cq}”=1  the  coefficients  that  must  be  selected  to  satisfy  inter¬ 
polation  conditions  or  constraints,  and  0  the  set  of  basis  functions.  The  interpolation 
constraints 


s(xi)  =  s,g  eR,  j  =  1  ,...,n  (6) 

are  enforced  on  the  linear  combination  (5)  to  ensure  that  the  approximation  matches 
exactly  at  locations  where  information  is  known.  The  radial  basis  functions  <f>{rj) 
depend  on  the  radial  parameter  r3  =  ||x  —  Xj||,  j  =  1,  ...,n,  which  implies  that  the 
value  of  the  radial  function  depends  only  on  its  distance  from  the  center  point  Xj 
-  and  is  thus  rotationally  invariant.  This  invariance  property  plays  a  key  role  in 
the  derivation  of  the  quadrature  weights  in  section  2.3.2.  The  choice  of  the  norm 
I'd  can  be  made  based  on  the  application,  but  is  most  often  the  Euclidean  two- 
norm.  The  definitions  of  the  radial  functions  differ  only  on  the  center  point,  so  the 
set  {0(| |x  —  Xj||)}”=1  is  a  collection  of  translations  of  the  same  function.  Using  the 
Euclidean  two-norm,  the  interpolant  takes  the  form 

s(x)  =  ^  Cj(ft  f  \J (ah1)  —  x^)2  +  (ad2)  —  x^)2  +  •  •  •  +  (ahn)  —  x^)2^  (7) 

j= i  ^  ' 
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Examples  of  RBFs  that  can  be  used  in  (7)  are  found  in  table  1.  In  determining  the 
coefficients  {cj}”=1,  the  interpolation  constraints  (6)  form  the  linear  system 

AC  =  G  (8) 


where 


0(||Xl  —  Xi||) 

0(1  |Xi  —  X2||)  • 

■'  0(1  |Xl  -Xn||) 

Cl 

A  = 

0(||X2  —  Xi||) 

0(1  |X2  —  X2||)  • 

'•  0(||X2-Xn||) 

,c  = 

C2 

_0(||xn-x1||) 

0(11  Xn  x2 1 1 )  • 

0(|  |xn  -Xn||)_ 

Cn 

and  G  = 


0X  l) 

0(x2) 


s(xjv) 


Thus,  A  is  an  n  x  n  matrix,  while  C  and  G  are  column  vectors  with  n  elements. 

For  many  choices  of  0(x),  the  matrix  A  can  be  shown  to  be  invertible  [7];  however, 
in  those  cases  where  A  is  singular,  (7)  can  be  regularized  with  polynomial  constraints: 


n  M 

s(x)  =  +  ^cfyp(x), 

j= i  i=i 


(9) 


cfBF7Ti(xj)  =  0,  for  l  =  1 M  (10) 

j=i 

to  guarantee  nonsingularity  [7].  That  is,  the  polynomial  terms  are  orthogonal  to  the 


19 


RBF  coefficients.  As  an  example,  for  x  G  M2,  M  —  (m  +  1  ){m  +  2)/2  and 

n(rj,  0  =  {1}  for  m  —  0 

ZfCffiO  =  for  m  =  1 

k(v,  0  =  {1,  V,  t  V2,  Vt,  £2}  for  m  =  2 

k(v,  0  =  {1,  V,  t  V2,  Vt,  £2,  V3,  V2t,  Vt2,  £3}  for  m  =  3 

7L(v,0  =  {hri,Z,V2,rt,e,-  ,riC~\C} 

The  inclusion  of  polynomial  terms  in  (9)  and  the  largest  degree  of  the  polynomial 
terms  included  (m)  can  have  impacts  on  the  accuracy  of  the  approximation  (see  for 
example  [10])  even  when  A  is  nonsingular.  Furthermore,  M  is  limited  by,  for  instance, 
(m+i)(m+2)  <  n  sjnce  invertibility  is  no  longer  possible  for  M  >  n  .  As  a  result  of  (9) 
and  (10),  the  linear  system  now  becomes 

AC  =  G  (11) 

where 

A  P  _  CRBF  _  G 

A  =  ,  C  =  ,  and  G  = 

PT  O(mxM)  Cp  O(mxi) 

Here,  A  and  F  are  the  same  matrix  and  vector,  respectively,  from  (8),  CRBF  and  Cp 
are  column  vectors  defined  by  CBBt  =  cRBF ,  j  =  1, ...,  N,  and  Cf  =  cf,  l  =  1, ...,  M, 
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as  from  (9),  while 


7Tm(x  i) 

ttm(x2) 

TM(xn) 

Thus,  P  is  an  n  x  M  matrix,  A  is  an  (n  +  M)  x  (n  +  M)  matrix,  and  C  and  G  are 
column  vectors  with  (n  +  M)  elements. 

Referring  back  to  table  1,  when  using  the  inhnitely  smooth  basis  functions,  the 
interpolation  matrix  A  is  positive-definite  (nonsingular).  Whereas,  the  monomial 
RBFs  (piecewise  smooth  functions)  applied  in  this  thesis  do  not  always  have  positive 
definiteness  [7].  As  a  result,  (9)  with  nonsingular  interpolation  matrix  A  is  used  to 
interpolate  g(r],£)  on  the  two-dimensional  coordinate  system  (77,  £)  of  section  2.2.2. 
Notice 

(m+l)(m+2) 

Ink{g)  ~  Ink{s)  =  c?BFlnk  (t  (y/ (v  -  Vj)2  +  (£  -  0)2)  )  +  0) 

(12) 


p  = 


7Tl(Xi)  7T2(Xi) 
7Ti(x2)  7T2(x2) 


TTilX,, 


7T2(X„ 


Dehne 


/  = 


jRBF 

P 


where  lfBt 


Ink  (0  (y/(v  ~  Vj)2  +  (Z  ~  Q2)) ,  3  =  1,2 ,...,n,  and  If  =  /(^(t?,  £)), 
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I  =  1,  2, M.  Since  AC  =  C  and  A  is  invertible, 

(a-'g)t  i 
=  gtw 

That  is,  W  =  I  is  the  solution  to  ATW  =  /.  Enumerate  the  elements  of  W 

as 

W  =  wBRF  wBBF  ■  ■  ■  wBBF  w\  w% 

Note  that  by  equation  (11)  only  the  first  n  elements  in  G  have  nonzero  values,  so 

n 

W»)  «  GTW  =  V  w^BFg(Vt,Q  =  ink(g)  (13) 

1=1 

2.3  Weight  calculations 

To  find  the  weights  {wBBF}r-=l  for  approximating  Ink(g),  the  linear  system  of 
(11)  must  be  constructed.  The  matrix  A  is  readily  populated  by  evaluating  each  of 
the  basis  functions  (RBFs  and  polynomials)  at  each  of  the  n  interpolation  points. 
The  construction  of  Ink(g),  on  the  other  hand,  requires  integrating  each  of  the  basis 
functions  (preferably  analytically)  over  the  planar  triangle  tabc- 

2.3.1  Planar  triangles. 

In  order  to  compute  the  weights  for  (13),  the  values  of  IBBF,  j  =  1,2,  ...,  n,  and 
If,  /  =  1,2,  ...,M,  from  (12)  are  required.  While  integrating  polynomial  terms  in  IF 
is  more  common,  this  section  discusses  a  way  to  obtain  the  values  IBBF  exactly. 
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2.3. 1.1  Integrating  Rotationally  Invariant  Functions  over  Right  Tri¬ 
angles. 

Let  O:  ujo  =  (vo,  io)  be  an  arbitrary  projection  point  in  the  (77,  £)  plane  defined  in 
section  2.2.2.  As  adapted  from  [2],  a  radial  basis  function  d>  ij]  —  r]0 )2  +  (£  —  £o)2) 
can  be  integrated  over  the  projected  triangle  tabc  £  hi k  with  vertices  A:  u>a  = 
(7 UiZa),  B:  ujb  =  (■ r)B,t,B ),  and  C:  urc  =  (rjc,£,c)  ( A ,  B,  C  being  projected  nodes) 
by  integrating  over  a  combination  of  six  different  right  triangles.  How  these  six 
triangles  relate  to  tabc  is  different  for  each  individual  choice  of  O.  Now  consider 
D :  <jl>d  =  (/7d,£d),  E :  ll>e  =  (r]F, £e),  and  F:  u>F  =  (r]F,iF)  to  be  the  orthogonal 
projection  of  O  onto  the  lines  through  sides  AB,  BC ,  and  AC  of  tabc,  respectively 
(see  figure  3).  Then  it  will  be  shown  that  the  integral  /(</>)  over  tabc  can  be  written 
as 

ItaboW)  '■=  ff  0  (■ V (v  -  Vo)2  +  (f  -  to)2)  didr]  =  soadUoad 0)  +  sodbUodb^) 

TABC 

+  SOBEhoBE (</0  +  SOEchoEci 0) 
+  SoEAhoFA {<t>)  +  SOCFhocF {<t>) 

That  is,  it  can  be  computed  as  a  linear  combination  of  the  integrals  over  the  right 
triangles  t0AD,  toDB,  t0BE,  toEC,  toFA,  and  tocF  with  corresponding  weights  s0ad, 
SODB,  sobe,  soec,  s0fa,  and  s0cf,  respectively,  equal  to  +1  or  -1. 

Consider  the  integral 


hoAD  (4>)  =  ff  0  ( V (v  -  Vo)2  +  (i  -  io)2)  didr]  (14) 

tOAD 

since  the  integrals  for  the  other  five  right  triangles  are  analogous.  Since  O  could  be 
located  anywhere  in  f it  will  generally  not  be  at  the  origin,  and  the  side  DO  will 


23 


Figure  3.  Example  divisions  of  a  planar  triangle  into  six  right  triangles,  (a)  Final 
orientation  of  an  individual  right  triangle  used  for  RBF  integration,  (b)  Example  division 
of  tabc  where  the  signs  soad,  sodb,  s0be,  s0ec ,  s0FA,  and  s0cf  in  equation  (19)  are 
all  positive,  (c)  Example  division  where  the  signs  so  AD ,  sofa,  and  socf  are  positive 
while  signs  sqdb,  sqbe ,  and  sqec  are  negative. 


not  align  with  the  ?/-coordinate  axis.  As  discussed  in  [37],  consider  moving  O  to  the 
(■ r /,  £)  origin  by  the  change  of  variables  (depcited  in  figure  4b) 


uj'  =  UJ  —  (jJq 


Next  consider  the  rotation  matrix 


cos# 

sin# 

and  R  1 = 

cos# 

—sin# 

R  = 

—sin# 

cos# 

sin# 

cos# 

where  #  is  the  counterclockwise  angle  between  the  positive  7/-axis  and  the  vector  in 
the  direction  of  the  side  DO.  Thus,  define  the  change  of  variables 

(j)  -  ?7o)cos0  +  (£  -  £o)sin6> 

~{V  -  770) sin#  +  (£  -  £o)cos# 

that  corresponds  to  translating  u0  to  (0,  0)  and  rotating  the  edge  DO  to  lie  on  the 
positive  77-axis  with  the  edge  AD  either  above  or  below,  but  perpendicular  to,  the 


uj"  =  Ruj'  =  R(u)  —  uj0)  = 
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//-axis  (see  figure  4c).  Note  that 


(n")2  +  (O2  =  tn  -  no)2  +  (Z-  to)2  (15) 


and 

cos9  —sin# 
sin#  cos  9 

=  cos2  9  +  sin2#  =  1 

So,  dr)d£  =  drj"d Therefore,  (14)  is  equivalently 

jjo  (V(#  -  #o)2  +  (£  “  £o)2)  d£dri  =  JJ  (j)  [y/Wf  +  (f")2)  drfd£" 

to  AD  toAD 

for  toAD  the  triangle  toAD  translated  and  rotated.  Now  with  O  at  the  origin  and 
DO  lying  on  the  positive  //-axis  with  AD  perpendicular  to  it,  rfn  =  rfA  and  t^b  =  0. 
Letting  old  =  \rfb\i  (15)  shows  that 

=  \nb\  =  \J (rfb)2  +  (&)2  =  V (VD  -  no)2  +  (Id  -  £o)2 

is  the  length  of  edge  DO.  Similarly  letting  (3oad  =  |£aI,  from  (15), 

bo  AD  =  \£a\  =  \/  ( nD  -  nA )2  +  (£a  -  £d)2 

is  the  length  of  edge  AD.  As  a  result,  the  integral  over  the  triangle  to  ad  is  now 

n  n  f*(A-D  r  jj/f 

JJ  4,  (vW+f?)  <T  =  J  J  °°  4>  (vW+FF)  W 

toAD 


dy" 

dz” 

dr i 

dr] 

dr\" 

dZ" 

&Z 

dZ 
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for  AD  above  the  //-axis  (see  figure  4d).  For  AD  below  the  //-axis  (figure  4e),  let 
rf"  =  fj"  and  where 


j\oAn  ,/(vW+CT)< 


dr,"  = 


(vVO2  +  (-r)2)  (--i)^'"////'" 
(\/ (v”')2  +  (CO2)  d?"dv"' 


In  any  case, 


p  p  ,  V  /*Oi£)  p  a  V 

Jj  (p  (V (h  -  ho)2  +  (£  -  £o)2)  dfdrz  =  J  j  HV V2  +  t2)d£dri  (16) 

tOAD 


Here,  ap  is  the  length  of  side  DO.  So,  and  op  would  be  the  lengths  of  sides 
EO  and  FO,  respectively.  Likewise,  Podb,  Poec ,  Pobe ,  Pofa ,  and  Pocf  could  be 
defined  in  a  similar  way  as  Poad- 


2. 3. 1.2  Some  Closed  Form  Integrals  of  RBFs. 

Note  that  many  commonly  considered  RBFs  can  be  integrated  exactly  over  right 
triangles  in  M2.  For  example,  given  a  right  triangle  with  base  a  and  height  P,  then 


Cr/ 


>0  Jo 


(r/2  +  £2)2 d^drj  =  a  (  3a4sinh  1  ( —  ]  +  P\/  a2  +  P2( 5a2  +  2 p2)  )  , 


a 


n—rj 

“  tf  +  eMdr, 


=  — —  a  ( 15a6sinh  1  (  —  )  +  pJ  a2  +  P2( 33a4  +  26a2 p2  +  8/34)  J  , 
336  \  \o . )  J 
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Figure  4.  (a)-(d)  is  the  progression  of  shifting  and  rotating  a  right  triangle  for  use  in  RBF 
integration  as  outlined  in  section  2. 3. 1.1.  Note  that  (a)  has  different  axes  for  scaling 


purposes.  In  the  case  of  u>£>  <  u>o  hi  (a),  (b)-(d)  would  look  upside-down,  resulting  in  (e). 
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and 


n—rj 

“  tf+eMdv 

a  ^105o;8sinh_1  (|)  +  /3 ^  oZ  +  /32(279o6  +  326a4/?2  +  200a2/?4  +  4806)) 
“  3456 


Each  of  these  are  defined  as  long  as  a  ^  0  (that  is,  when  the  vertices  of  the  right 
triangle  are  not  all  three  collinear)  [2],  Hence,  these  closed  form  integrals  of  monomial 
RBFs  can  be  used  to  evaluate  such  integrals  as  (16). 


2.3. 1.3  Integration  over  Arbitrary  Triangles  via  Integration  over 
Right  Triangles. 

To  integrate  radial  basis  functions  over  planar  triangles,  first  consider  a  right 
triangle  txYZ  hi  the  (77,  £)  plane  made  up  of  the  line  segments  from  X\  (rjx,£x)  to  Y, 
Y:  ( r}y,£y )  to  Z,  and  Z:  (r)x,£z)  to  X.  Suppose  that  (f>(r),£)  is  integrable  with  respect 
to  both  77  and  £  and  define 

Mv,0  =  +\f  <t>(v,Z)dx  and  Mv,0  =  ~  f 


(this  holds  for  all  the  radial  basis  functions  used  this  thesis,  as  well  as  multivariate 
polynomials).  Then, 

d  d 

nv,0  = 

Dehne 


/ 

1 

V 

Mv,0 

V 

J 

Mv,0 
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so  that  the  line  integral  of  <f>  along  the  line  from  X  to  Y  can  be  written  as 


LIxy{$)  '■= 


1 

/ 

\ 

/  $ 

Vx 

+  t 

Vy 

-  Vx 

Vy 

-  Vx 

J 

0 

V 

1 

Zy 

1 

Uy 

i 

) 

Zy 

- 1 

Uy 

1 

i 

/ 

\ 

/  $ 

Vy 

+  t 

Vx 

-Vy 

Vx 

-Vy 

J 

0 

V 

1 

Uy  j 

Zx 

1 

Uy 

1 

J 

Zx 

- 1 

Uy 

1 

~ LIyX{ $) 


Since  the  boundary  of  txYZ  is  a  closed  and  piecewise  smooth  contour,  and  so  long  as 
X,  Y,  and  Z  are  oriented  counterclockwise  in  the  plane,  Green’s  theorem  [29]  allows 
the  integral  of  txYZ  to  be  written  as 


ItxYZ  (0)  =  LIXy{$)  +  LIyz{®)  +  LIZX{$) 


Note  that  if  < £)  =  1,  then 

Area  (tXYZ)  =  hXYZ(0) 

^txYZ  (1)  >0 

Also  note  that  -§^{\v)  —  §0  =  1-  Continuing  with  a  counterclockwise  orientation, 


4xxz(l)  - 

( 

1 

Uy 

rH|W 

1 

1 _ 

\ 

+  -h-fyz 

/ 

1 

Uy 

'-HCN 

1 

1 _ 

\ 

+  hUx 

/ 

1 

Uy 

rHlCM 

1 

1 _ 

\ 

1_, 

l„ 

V 

U 

/ 

V 

toi 

/ 

V 

U 

/ 

1 

1 

Uy 

1 

1  ^ 

Vz  -  Vx 

2 

1 

* 

c- 

1 

p-  , 

- 1 

><: 

Uy 

1 

Uy  ( 
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With  a  clockwise  orientation,  the  direction  (and  therefore  sign)  of  the  path  changes: 


^xyz(^)  LIXZ 


/ 

1 

rH|CM 

1 

\ 

V 

. 

/ 

/ 

1 

■Uy> 

i-HlCN 

I 

\ 

( 

1 

■Uy> 

rHlCM 

1 

\ 

+  hT^y 

+  LIyX 

i„ 

1_, 

V 

i 

ICM 

) 

V 

27/ 

/ 

1 

1 

■U j 

1 

1  ^ 

Vz  -  Vx 

2 

1 

S' 

1 

S'  , 

iz  -  ix 

Therefore,  let 


Sxyz  :=  sign 


1 

>1 
■U j 

1 

1  ^ 

1 

>< 

S' 

I 

N 

S' 

\ 

1 - 

1 

* 

1 

1 

(S) 

Uyi 

/ 

(17) 


where  sxyz  =  1  if  the  vertices  of  txYZ  are  oriented  counterclockwise,  and  sxyz  = 
— 1  if  they  are  oriented  clockwise.  It  is  apparent  that  if  the  vertices  are  oriented 
counterclockwise,  then  sxyz  =  1  and 


SXYzItXYZ{4> )  =  I^XYZ  (<t>)  =  LIxy{$)  +  LIyz{$)  +  Llzx^f) 


Whereas  vertices  oriented  clockwise  yields  sXyz  =  —  1  and 


SXYzItXYZ{4> )  -  (-1  )hXYz{4>)  -  (-l)(^vz(<f>)  +  T/Zy($)  +  LIYX{$)) 

=  (-1  ){-LIzx{$)  ~  LIyz ($)  -  LIXY{<$>)) 
=  LJzxi^)  +  +  h/yy  (T) 


In  either  case, 


sxYzhxvz (0)  —  +  T/yZ($)  +  L/Zy(<h)  (18) 
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Recall  that  the  task  at  hand  is  to  show  that  the  integral  over  an  arbitrary  right 
triangle  can  be  evaluated  by  computing  integrals  over  six  right  triangles  with  common 
vertex  ().  Parts  of  these  right  triangles  may  reside  within  tabc,  while  other  parts 
may  reside  outside.  A  combination  of  the  six  right  triangles  should  create  the  area 
for  the  ABC  triangle.  By  adding  or  subtracting  the  integrals  of  the  six  right  triangles 
in  question,  it  is  claimed  that  the  integral  value  IfABC{4>)  over  Tabc  can  be  found  by 
evaluating 

IfABcirf*)  ^ABc{^OAdC-oAd{C ft)  A  SODBhoDB^)  A  SOBElt0BE  {4>)  A  SOEcItoEC  (</>)+ 

SOCfUocf^)  +  sOFaA0fa  (</*))  (19) 

for  signs  s0ad ,  s0db,  s0be,  s0ec,  s0fa,  and  s0cf  defined  as  in  (17)  by  replacing 
XY Z  and  preserving  vertex  order.  Depending  on  the  orientation  of  the  six  right 
triangles  and  tabc ,  the  signs  may  all  be  1,  or  only  some  may  be  equal  to  1.  As  seen 
in  figure  3,  (b)  depicts  a  tabc  such  that  all  the  signs  are  positive,  while  (c)  depicts  a 
tabc  such  that  the  signs  soad,  sofa,  and  socf  are  positive  while  signs  sodb,  sobe, 
and  soec  are  negative.  Regardless,  the  orientation  of  the  vertices  A,  B ,  C,  D,  E,  F, 
and  O  is  the  deciding  factor  on  delegating  the  values  for  the  signs. 

Recall  that  Green’s  theorem  [29]  is  based  on  a  boundary  (path)  integral  and  that 
LIqa  =  —LI Ao-  With  all  this  information,  the  following  works  backwards  to  achieve 
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(19)  [38].  Similar  to  figure  3(c), 


SABc{sOADhOAD  {4>)  +  SODBhooB  (0)  +  sOBEhOBE  (0)  + 

SoEchosc ($)  +  sOCf4ocf(</>)  +  sOFAltoFA (0)) 

=  sabc(LIoa(*&^  +  LIad{$)  +  £/,do(<£)+ 

LIod{^)  +  LIdb{Q)  +  ^/_bo(®)  + 

LIob{^)  +  LIbe{^)  +  ^^eo(<^>)  + 

L/oe($)  +  L/sc($)  +  -^^co(®)+ 

L/oc($)  +  LIcf{Q)  +  ^^fo(<^>)  + 

L/of($)  +  £/fa($)  +  ^o($)) 

=  sabc^-Tad^)  +  LIdb{$)  +  ^/bb(<1>)  + 

L/sc($)  +  L/cf($)  +  L/fa($)) 

=  Sabc{LIab{$)  +  -^-TbC^)  +  -^^Ca(^))  =  If abc  (0) 

Now  that  the  integrals  of  RBFs  for  individual  triangles  can  be  evaluated  exactly,  the 
following  puts  all  the  previous  Chapter  If  information  together  to  yield  the  desired 
quadrature  weights  for  Tg(f). 

2.3.2  Converting  Weights  from  Planar  to  Spherical. 

Recall  from  sections  2.2.1  and  2.2.2  that  the  nearest  neighboring  nodes  for  each 
triangle  rfc  G  T  (that  is,  the  points  in  J\fn)  are  projected  onto  the  (77,  £)  plane,  and 
then  a  transformation  allows  the  projected  nodes  to  be  treated  as  though  they  exist 
in  two-dimensional  space.  With  intermediate  quadrature  weights  from  (13)  in  hand, 
the  task  of  relating  the  weights  computed  for  each  tangent  plane  to  the  ones  required 
for  computing  surface  integrals  on  the  sphere  remains. 
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Considering  the  change  of  variables  (3)  and  (4)  from  section  2.2.2,  inverse  trans¬ 


formations  mapping  points  in  0^  back  to  Snk  can  be  realized  by 


x  = 


y  = 


z  = 


VP 

\J  p1  +  rf  +  £2 

jp 

\J  P1  +  T]2  +  £2 


(20) 


in  the  case  of  x2M  +  y2M  =  0.  Similarly,  when  the  triangle  midpoint  Xm  is  not  located 
on  the  z-axis  (x2M  +  y2M  ^  0),  the  inverse  transformation  is 


xM  (rjzM  +  Py/xM  +  vh)  ~  p£vm 


x  = 


y  = 


y/p2  +  v2  +  i2VS:M  +  Vm 

Hm  (vzm  +  Py/xh  +  yh)  +  pixM 


(21) 


z  = 


y/p2  +  V2  +  i2y/x2M  +  y2M 

~Vy/x2M  +  y\i  +  p~zm 

y/ p2  +  p2  +  £2 


That  is,  Tfc  G  T  is  parameterized  locally  by  the  parameters  p  and  £.  Given  a  surface 
parameterization  x(rj,  £),  y(rj,  £),  z(jj,  £),  it  is  true  by  definition  [39]  that 


f(x,y,z)dS=  //  f{x(r],t),y{r],t),z{r],t)) 


|-x(r/,0  x  j^x(p,£) 


dyd£ 


Tk 


Tk 


In  either  case  of  (20)  or  (21),  and  after  much  simplification, 


|-X(/7,0  x  |fx(tpO 


P 


2  (p2  +  P2  +  £2)  2 
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So, 


f(x,  y,  z)dS  —  //  f(x(r],Z),y(ri,£),z(ri,t)) 


P 


Tk 


Tk 


C p 2  +  v2  +  £2) 


■drjd £ 


Defining 

P3 

g(v,0  =  f(x{v,£),y{v,£),z{v,0)  — — - — — i, 

(, p 2  +  t  +  i )2 

the  integral  over  rk  can  be  approximated  as  discussed  in  the  preceding  sections  to 
give 


f(x,  y,  z)dS  «  Y 

3= i 


wfBFf(x(r),  C),y(v,  f)>  z(v,  0) 


(p2  +  rf  +  £2)  2 


(22) 


Notice  that  the  parameterization  for  each  triangle  rk  is  different,  so  the  respective 
parameterization  will  now  be  indexed  as  r)k  and  £*..  Considering  the  surface  integral 
over  all  of  the  sphere  S, 


I< 


Mf)  =  £  f(Vk,£k) 


k= i 


p 


(p2+vt+ek)^ 


drikd^k 


I< 

££< 

fc=i  j=i 


^BF)3f((Vkb 


P 


(p2  +  (%)?  +  &)?)* 


Here,  {(wFRF).,}’-=1  is  the  set  of  weights  in  (22)  for  triangle  rk  and  {(%).?'}  j=u 
{(£fc)j}j=i  represent  the  n  nearest  neighbors  in  Sn  to  Xm  after  the  parameteriza¬ 
tion  is  applied. 

Let  K-i ,  i  =  1,2,...,  At,  be  the  set  of  all  pairs  (/c,j)  such  that  {(Vk)j,  (£k)j)  ^ 
(. Xi,yi,Zi ).  That  is,  each  node  (. x^y^Zi )  maps  to  at  least  one  respective  (r}k,£k) 
plane.  Indeed,  this  easily  occurs  more  than  once  as  each  triangle  rk  projects  n  nearest 
neighbors.  The  larger  n  is,  the  more  times  each  node  will  be  projected  into  separate 
tangent  planes.  So  each  node  ( Xi,yi,Zi )  has  quadrature  weights  from  a  number  of 
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triangles,  and  the  sum  of  those  weights  from  all  the  planar  regions  it  was  projected 
into  is  called  Wt .  Then  the  integral  Ts(f)  can  be  written  as 


N 

Ml)  ~  £ 

1=1 

N 


£ 

Ak,j)eKi 


KBFh 


P 


(P2  +  fafc)?  +  &)?) 


=  '52wif(xi,yi,zi)=is{f) 


f(x i,  Vi,  -  ) 


(23) 


As  with  previous  parts,  the  above  calculations  are  referenced  from  [2],  Equation  (23) 
is  the  cumulation  of  the  works  of  Chapter  II,  and  is  the  numerical  approximation  of 
the  integral  (i.e.  quadrature)  of  some  function  f(x,y,z)  over  the  surface  of  a  sphere. 

The  accuracy  of  (23)  will  depend  on  the  accuracy  of  the  RBF  interpolation  scheme 
employed  for  each  planar  region  Q*,  [2] .  Over  scattered  data  sets,  the  accuracy  is  often 
reported  relative  to  the  Ell  distance  [40]  (a.k.a.  mesh-size  or  mesh  norm  [25]).  The 
fill  distance  is  defined  based  on  the  situation.  In  some  cases,  the  mesh  size  is  the 
(largest)  distance  between  any  two  points.  In  the  case  of  the  sphere  used  in  this 
chapter,  the  fill  distance  can  be  described  as  the  diameter  (in  arc-length)  of  the 
largest  empty  spherical  cap  between  two  points  on  the  sphere  surface.  Whereas  on 
the  tangent  planes,  the  fill  distance  is  the  diameter  of  the  largest  empty  planar  circle 
in  the  convex  hull  of  the  point  set.  In  other  words, 


h\rnk  ■■=  2sup:rgQfe dist (x,  M) 

where  J\f  is  the  set  of  points  in  a  bounded  domain  and  hj^^k  is  the  diameter 
(not  radius)  of  the  largest  empty  circle  between  any  two  points  in  Q&  (adapted  from 
[2,  25]). 

When  projecting  the  points  on  the  sphere  to  an  (r/^,  £fc)  plane,  the  sphere  fill 
distance  may  become  distorted  as  per  (20)  and  (3).  As  this  affects  the  accuracy  of 
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the  approximation  integral  from  (23),  verifying  that  the  distortion  will  not  decrease 
the  accuracy  significantly  for  a  small  change  in  sphere  fill  distance  is  important  [2], 

2.3.3  Some  Test  Cases  Over  a  Sphere  Used  in  [2]. 

A  total  of  six  test  cases  are  presented  in  [2],  of  which  two  in  particular  are  re¬ 
displayed  here  for  the  purpose  of  comparison  in  section  4.2.  For  the  unit  sphere  §>2, 
these  are 


fi(x,  y,  z)  =  1  +  x  +  y2  +  x2y  +  xA  +  y5  +  x2y2z 2 

~  .  .  1  +  sign (—9a;  -9 y  +  9z) 

J2{x,  y,  z)  = - - - 


whose  exact  surface  integrals  are 


X§2(/i) 

2s*  (/2) 


216vr 
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Consider  the  application  of  the  method  discussed  in  this  chapter  with  0(r)  =  r  , 
m  —  7  (degree  of  polynomial),  and  n  =  80  (nearest  neighboring  nodes).  Note  that 
fi(x,y,z)  G  (^“(S2)  with  convergence  rates  from  the  RBF  quadrature  method  dis¬ 
played  in  figure  5.  Meanwhile,  f2(x,y,z)  is  a  discontinuous  yet  bounded  function 
whose  convergence  rates  from  the  RBF  quadrature  method  are  displayed  in  figure 
6.  The  relative  errors  reported  in  figures  5  and  6  are  the  maximum  error  over  1000 
random  rotations  of  each  node  set.  The  minimum  energy  (H)  nodes  were  generated 
by  [2]  using  the  methods  from  [32],  which  appear  in  [25].  The  minimum  energy  (W) 
nodes  were  generated  from  [5].  It  was  shown  in  [2]  that  various  test  integrands  in 
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Error  in  Quadrature 
Over  the  Sphere  Surface 
Test  Function  ,/i 
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Figure  5.  Relative  error  in  RBF  quadrature  for  test  function  fi(x,y,z)  over  M2  as 

in  [2]. 
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Figure  6.  Relative  error  in  RBF  quadrature  for  test  function  f2(x,y,z)  over 

in  [2]. 


as 


applied 


applied 
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C°° (§2)  resulted  in  a  convergence  rate  of  0(N~3'5),  as  illustrated  in  figure  5.  Since 
f2(x,y,z)  is  a  discontinuous  function,  its  expected  convergence  rate  should  be  worse 
than  that  of  a  smooth  function. 

ft  was  noted  in  section  2.2.3  that  as  the  degree  of  interpolating  polynomial  m 
increases,  then  the  accuracy  of  the  interpolation  should  conceptually  increase  as  well. 
An  application  of  varying  m  is  depicted  in  figure  7  verifying  this  concept.  However, 
just  as  the  accuracy  increases  with  increasing  m  and  n,  so  too  does  the  CPU  time  for 
computing  the  quadrature  weights  increase.  As  a  result,  balancing  between  accuracy 
and  computation  cost  is  dependent  on  the  application. 

For  fixed  m  —  7,  n  —  80,  and  <f>(r)  =  r7,  the  total  cost  in  computing  the  RBF 
quadrature  weights  over  §2  is  illustrated  in  figure  8  [41].  Not  only  is  the  spherical 
quadrature  method  in  [2]  accurate  at  0(iV~3'5)  (for  smooth  integrands  /),  the  method 
is  computationally  inexpensive  at  O(NlogN)  operations  and  O(N)  memory  usage. 
When  expanding  the  computation  of  weights  beyond  one  core,  parallel  scalability 
with  number  of  cores  was  also  observed  (see  figure  8c). 

With  the  conclusion  of  Chapter  II,  the  following  chapter  adapts  the  numerical 
quadrature  method  for  a  sphere  surface  to  a  quadrature  method  for  a  smooth  closed 
surface.  Because  much  of  the  methodology  remains  the  same  in  approximating  a 
smooth  versus  sphere  surface,  a  separate  chapter  was  dedicated  to  spherical  quadra¬ 
ture  as  a  build-up  to  the  more  recent  research  on  surface  quadrature. 
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Log  Base  10  of  Error 
in  the  Surface  Integral  of 
f(x,y,z)  =  cos  (fz) 
with  N  =  22500  (RBF  r‘) 


CPU  Time  (Seconds) 
N  =  22500  (RBF  r7) 


Figure  7.  (Adapted  from  [2])  For  some  example  function  f(x,y,z)  =  cos(^z)  with  a  given 
number  of  nodes  and  radial  basis  function  <fi(r)  =  r7,  (left)  indicates  the  quadrature  error 
with  respect  to  the  number  of  nearest  neighbors  n  and  degree  of  polynomial  m  used  in 
equation  (9),  (right)  indicates  the  CPU  time  (in  seconds)  for  computing  the  quadrature 
weights  with  respect  to  n  and  m.  The  black  dashed  line  indicates  the  constraint 
^m+1^m+2^  <  rb  such  that  the  interpolation  problem  becomes  singular  when  the  constraint 
is  not  met.  These  plots  were  generated  in  machines  with  dual  Intel  Xeon  E5-2687W  3.1 

GHz,  8-core  processors. 
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Figure  8.  (Adapted  from  [2])  Timing  results  for  computation  of  the  RBF  quadrature 
weights  for  the  surface  area  of  M2  as  applied  in  [41],  The  spherical  quadrature  method  has 
rate  0(N )  for  computation  cost  and  memory  usage  for  at  least  millions  of  nodes,  while  a 
cost  of  O(NlogN)  operations  is  expected  for  when  the  Delaunay  triangulation  and  nearest 
neighbor  search  dominate  the  computation  [41],  Parallel  scalability  with  number  of  cores 
was  also  observed.  These  plots  were  generated  in  machines  with  dual  Intel  Xeon 
E5-2687W  3.1  GHz,  8-core  processors. 
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III.  Quadrature  for  Smooth  Closed  Surfaces 


As  the  previous  chapter  focused  on  approximating  the  integral  of  some  scalar  func¬ 
tion  f(x,  y,  z )  over  the  surface  of  a  sphere,  this  chapter  too  approximates  the  integral 
of  some  scalar  function  f(x,y,z),  but  over  a  generalized  closed  smooth  surface.  The 
surface  S  can  either  be  predefined  implicitly  by  h(x,  y ,  z)  =  0,  predefined  explicitly  by 
a  parameterization  of  coordinates  x(u,  v),  y(u,  v),  and  z(u,  v ),  or  simply  be  expressed 
through  a  set  Sn  of  N  points  in  M3.  However,  if  the  surface  is  not  predefined  by 
an  equation  and  is  only  described  by  a  set  of  quadrature  nodes  Sn,  then  the  surface 
shape  appears  as  a  clustering  of  N  points.  This  being  the  case,  the  remainder  of  this 
paper  assumes  that  in  the  case  of  only  Sn  being  given  (no  surface  equations),  then  a 
triangulation  T  with  triangles  {rfc}f=1  is  also  given  whereby  the  nodes  Sn  make  up 
the  vertices  of  the  triangles  -  just  as  in  Chapter  II. 

A  process  for  discretizing  a  surface  described  by  h(x,  y,  z)  —  0  using  N  nodes  and 
an  associated  triangulation  is  described  in  the  section  3.1.  While  a  surface  described 
explicitly  by  a  parameterization  of  coordinates  x(u,v),  y(u,v),  and  z(u,v )  can  be 
expressed  by  discrete  nodes  and  a  triangulation,  this  thesis  does  not  go  into  details 
about  how  this  is  done.  This  is  a  topic  left  for  future  consideration. 

Just  as  with  Chapter  II,  this  chapter  seeks  to  approximate  the  integral  of  f(x,  y,  z), 
but  over  some  smooth  closed  surface  S  as  in  (1): 


Uf)  ■■= 


N 

f(x,  y,  z)dS  «  ^2  wif(xii  y»  zi ) 

i= 1 


with  quadrature  weights  WW:  {Wi\S=l  for  N  nodes.  The  organization  of  this  chapter 
is  similar  to  that  of  the  previous  chapter,  in  that  it  will  walk  through  the  steps 
for  determining  the  quadrature  weights  in  (1).  As  various  steps  are  similar,  many 
references  to  Chapter  II  will  be  made  while  the  differences  are  expanded  upon.  The 
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Figure  9.  Example  surface  of  revolution  of  a  Cassini  oval  which  has  been  discretized  using 
distmeshsurface.  The  quadrature  method  of  Chapter  III  focuses  on  approximating  the 

integral  over  each  individual  triangle. 

material  for  this  chapter  was  derived  from  [42] . 

3.1  Triangulation 

Consider  a  smooth  closed  surface  S  defined  by  h(x,  y,  z )  =  0.  While  the  sur¬ 
face  is  defined,  the  quadrature  method  for  this  chapter  requires  quadrature  nodes 
{(xi,  yi,  Zi)}^=l  e  Sn  and  corresponding  triangulation  T.  Suppose,  however,  an  ap¬ 
plication  requires  both  Sn  and  T  to  be  composed  (only  given  h(x,y,z)  =  0).  This 
section  focuses  on  a  MATLAB  program  called  distmeshsurface  by  Persson  and 
Strang  [43],  which  does  just  that. 

Distmeshsurface  takes  a  surface  function  h(x,  y,  z)  =  0  and  outputs  equilibrium 
nodes  based  on  its  algorithm,  as  well  as  the  necessary  triangulation  T.  Recall  the 
earlier  description  for  fill  distance  in  section  2.3.2.  The  mesh  size  for  this  algorithm 
can  be  described  as  the  closest  distance  between  any  two  points  (Euclidean  two  norm). 
While  the  full  algorithm  can  be  found  at  [43],  a  summary  of  its  mechanism  can  be 
described  by  the  following: 

1.  Fill  a  user-defined  domain  with  points  separated  by  a  given  mesh  size,  and 
connect  the  uniformly  distributed  points  with  edges  to  create  a  tessellation. 
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2.  Points  whose  edges  are  completely  outside  the  surface  S  (edges  not  crossing 
through  S)  are  discarded,  and  the  new  grid  is  re-triangulated. 

3.  Following  the  re-triangulation,  the  points  are  moved  around  (while  ensuring  they 
remain  approximately  on  S)  to  make  the  edge  lengths  as  uniform  as  possible 
through  a  repelling  force  dependent  on  the  current  edge  lengths. 

4.  Steps  2  and  3  repeat  until  the  edge  lengths  are  uniform  within  a  set  tolerance. 

The  algorithm  then  outputs  quadrature  nodes  6  I3  and  triangulation  T,  as  in 
figure  9. 

3.2  Locating  Projection  Origin(s) 

Now  with  quadrature  nodes  Sn  and  triangulation  T,  the  next  step  in  this  thesis’ 
quadrature  method  is  to  reduce  the  problem  dimensionality  from  three  to  two.  As  in 
Chapter  II,  the  integral  over  each  surface  triangle  is  considered  separately.  However, 
the  flat  triangles  in  T  are  only  approximations  to  the  surface,  not  actually  part  of  the 
surface  whose  edges  are  geodesics.  As  a  result,  the  projection  and  change  of  variables 
for  each  triangle  {t/c}(L1  e  T  must  be  made  so  that  the  sum  of  the  projected  areas  do 
not  exceed  the  surface  area  of  S.  This  is  done  by  first  locating  a  projection  “origin” 
x<3,  or  point  through  which  a  portion  Snk  of  the  surface  is  projected  onto  the  plane 
through  the  triangles  vertices  (see  figure  12). 

To  start,  consider  the  three  edges  of  triangle  Tabc  £  T:  AB,  AC,  and  BC  with 
vertices  x^:  (xa,Ua,  za),  xs:  (xb,  Ub,  zb),  and  xc:  ( xc,yc,zc )•  Now  consider  an 
adjacent  triangle  tabe  with  third  vertex  xg:  (xB,  He,  ze)  (see  figure  10).  Triangles 
tabc  and  tabe  share  the  edge  AB.  Also  note  that  the  planes  containing  these  tri¬ 
angles  each  have  their  own  respective  normals  to.abc  and  tciabe  of  tabc  and  tabe, 
respectively. 
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Figure  10.  Taking  a  region  Snk  about  tabc ■>  this  illustrates  the  computation  of  the  edge 

normals  for  edges  AB1  CA,  and  BC. 


Let 


nAB  =  -  (nABC  +  Agn(nTABCnABE)nABE) 

describe  the  average  of  the  two  normals  that  lies  on  the  edge  AB.  The  same  can  be 
done  similarly  to  define  n Ac  and  n bc,  as  in  figure  10. 

Define  the  “cutting  plane”  through  the  edge  AB  to  be  the  plane  containing  the 
edge  AB  such  that  the  plane  is  parallel  to  nAs-  Consider  the  cutting  planes  passing 
through  the  edges  AB  and  AC  of  triangle  tAbc ■  Also  consider  the  vectors  wAB  (which 
points  from  to  x^)  and  Vqa  (which  points  from  x,^  to  x^).  Note  that  the  vectors 
nAB  and  vAB  are  both  parallel  to  the  cutting  plane  through  the  edge  AB.  Likewise, 
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(a)  (b) 

Figure  11.  (a)  Taking  the  cross  product  of  n ab  and  vab  yields  the  perpendicular  vector 
nOAB-  (b)  Taking  the  cross  product  of  n oab  and  nocA  yields  the  perpendicular  vector 

V(M- 


the  vectors  n ca  and  vCa  are  both  parallel  to  the  cutting  plane  through  the  edge  AC . 
In  order  to  mathematically  define  the  cutting  plane  through  the  edge  AB ,  a  normal 
vector  to  the  plane  is  required.  From  the  vector  cross  product, 


noAB  =  n ab  x  vab- 


where  the  subscript  O  indicates  that  the  projection  origin  xo  also  falls  in  the  AB 
cutting  plane  (see  figure  11a).  Similarly, 


nOCA  —  n-cA  x  vCA 


nOBC  —  nBC  X  WBC 


From  n oab  and  n oca,  a  vector  voa  (pointing  in  the  same  direction  as  the  vector 
from  xq  to  x^)  can  be  computed  via 


voa  —  n0As  x  n OCA 
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as  shown  in  figure  lib.  The  line  in  the  direction  of  vpA  and  passing  through  can 
be  parameterized  by 


x0  =  xA  +  tw0A  (24) 

such  that  x<3  is  the  point  along  this  line  that  intersects  the  cutting  plane  for  the  edge 
BC.  This  cutting  plane  can  be  represented  by 


n obc  ■  (x  -  xB)  =  0 


(25) 


where  •  is  the  vector  dot  product.  Solving  for  t  from  (24)  and  (25)  yields 

t  _  nOBc  •  (xg  -  xA) 

Hose  ■  vOA 


Hence,  the  projection  origin  used  for  triangle  tabc  is 

,  n obc  ■  (xs  -  xA) 

x0  =  xA  H - V0A 

nOBC  '  vOA 

While  the  above  can  certainly  be  used  as  the  projection  origin  for  projecting  the 
triangle  area  into  a  plane,  notice  that  the  node  xp  was  not  included.  If  a  more 
symmetric  version  of  the  projection  origin  (with  regard  to  all  three  vertices  of  tAbc ) 
is  preferred,  it  can  be  written  as 


xo  —  xM  +  —  [(n obc  •  (xb  —  xA))v<x4  +  (n oca  ■  (xc  —  xb))^ob  +  (n oab  ■  (xA  —  xp))v0o] 

(27) 
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where 


V  =  ((xc  -  xB)  x  nBC)  ■  [((xB  -  xA)  x  nAB)  x  ((xA  -  xc)  x  nCA)] 

=  ((xA  -  xc)  x  nCA)  •  [((xc  -  xB)  x  nBC)  x  ((xB  -  xA)  x  nAB)] 

=  ((xB  -  xA)  x  nAB)  •  [( (xA  -  xc)  x  xiCa)  x  ((xc  -  xB)  x  nBC)] 

3.3  Projecting  the  Nearest  Neighboring  Nodes 

Just  as  with  section  2.2.1,  the  n  nearest  nodes  A/"n:  {xj}"=1  G  Sn  to  each  triangle 
are  needed  for  numerical  approximation  via  RBF  interpolation.  When  projecting 
these  neighbors,  care  must  be  taken  to  avoid  projecting  points  for  which  the  vector 
from  the  projection  origin  Xo  to  the  point  being  projected  (x  G  A fn)  is  not  an  angle 
more  than  90°  from  the  normal  n asc-  Similar  to  section  2.2.1,  points  at  this  angle 
result  in  singularities  in  the  following  mathematics.  Therefore,  points  having  these 
conditions  are  not  utilized  for  projection.  That  is,  consider  the  projection  of  a  point 
x  G  Sak  to  the  plane 


nABc  •  (x  -  x.u)  =  0 

containing  triangle  tABc,  for  xm  the  midpoint  of  tABc ■  The  projection  illustrated  in 
figure  12  can  be  performed  by  finding  the  intersection  of  the  line 

X  =  Xo  +  t(x  -  Xo) 


47 


where  Xo  is  taken  from  (27)  and  x  is  the  projection  of  x  e  Sn  onto  the  Tabc  plane. 
Plugging  the  x  equation  into  the  equation  for  the  plane  yields 

^  _  H-abc  '  (xm  —  xp) 
n abc  ■  (x  -  x0) 


so  that 


x  =  x0  + 
=  x0  + 


n  abc  ■  (xm  -  xo),  x 

- - - -  x  -  x0 

n abc  •  (x  -  x0) 

nABc(xM  -  xQ)  +  n%c(yM  -  j/g)  +  n(ABC(zM  -  zp)  / 

nABc(x  -  xo )  +  «abc(2/  -  2/0)  +  -  *o) 


Xo) 


(28) 


wherenASC:  (n$£c,  n%c,  n%c),  xM:  (zM,  yM,  zM),  x:  (z,  y,  2),  and  xG:  (z0,  2/0,  ^o). 
Another  example  of  this  process  is  illustrated  by  the  points  in  figure  13a  being  pro¬ 
jected  as  in  figure  13b.  Note  that  (28)  differs  from  the  spherical  case  in  section  2.2.1 
only  in  that  Uabc  is  no  longer  the  vector  %  through  the  origin,  while  xq  is  no  longer 
at  the  origin.  Following  the  same  procedure  as  with  Chapter  II,  the  projected  nodes 
{Xj}”=1  will  then  be  used  to  define  a  two-dimensional  coordinate  system  on  the  Tabc 
plane. 


3.4  Defining  a  Two-Dimensional  Coordinate  System 

Recall  from  section  2.3.2  that  in  order  to  simplify  the  transformation  from  three- 
dimensional  coordinates  to  two-dimensional,  points  in  the  three-dimensional  coordi¬ 
nate  system  were  multiplied  by  a  rotation  matrix 
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Figure  12.  Projection  of  region  Sak  onto  planar  region  by  projection  origin  xq. 
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where  xM:  (xm,Vm,  %m)  was  the  spherical  triangle  midpoint,  p  =  sj x2M  +  y2M  +  z2M 
and  7  =  \Jx2M  +  y2M.  This  rotation  matrix  rotated  any  triangle  midpoint  Xm  to  the 
z-axis  at  (0,0,  p)  while  also  preserving  length  and  area.  For  an  arbitrary  smooth 
closed  surface  S,  the  transformation  is  more  complex  as  the  projection  origin  X£  for 
each  77  is  not  necessarily  located  at  the  origin  (0,  0,  0). 

In  order  to  achieve  a  similar  result  as  was  done  with  the  spherical  quadrature  in 
section  2.3.2,  the  three-dimensional  coordinate  system  is  translated  such  that  xo  is 
moved  to  the  origin.  Thus,  define  the  coordinates 

/// 

X  =  X  —  Xo- 

for  x  defined  in  (28).  While  the  surface  S  may  not  be  symmetric  about  the  origin 
as  with  the  sphere,  the  points  xw  can  still  be  rotated  about  the  origin  to  align  n abc 
with  the  vertical  axis.  As  a  result,  the  rotation  matrix  becomes 


M\2  — 


nABCnABC 

nABCnABC 

~7n 

'Jnpn 

' "Ynpn 

pn 

Ay) 

„(*) 

nABC 

nABC 

o 

In 

In 

(x) 

( y ) 

(z) 

n  ABC 

n  ABC 

nABC 

Pn 

Pn 

Pn 

where 


In  = 


2 


Pn  — 


2 


for  n abc'-  (■ nABC->  nABC->  nABC )•  Notice  that  is  of  the  same  form  as  rotation  matrix 
M,  only  differing  by  the  same  reasons  that  the  projection  into  the  plane  differs,  as  in 
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section  3.3.  Rotating  the  translated  coordinate  system, 


x"  =  M12x"' 


=  Mi2(x  -  x0) 


n-ABC  ■  (XM  —  xo)  ,,r  / 
—  - 7 - —  Mi2(^X 

nABC  ■  (x  -  Xo) 


Xo) 


tlABC  •  (XM  —  Xp) 

IlABC  •  (x  -  Xo) 


(x) 

nABCnABC 


ABc(X  Xo)+nABC’*ABG 


_ _ _ nABc(y  Vo)  (('nABc')~  +  l'nABc'l~')(z  z°'> 

(nABc)2  +  (nABc)2\/(.n%c)2  +  (nABc)2  +  (nABc)2 
-nABc(x-xo)+n{ABc(y-yo) 

\/  (nASc)2  +  (nABc)2 

r  ^  71  ~  ^  71  "  /77  TT 


constructs  a  new  three-dimensional  coordinate  system  such  that  Xm,  x,  and  xo  are 
the  same  as  in  section  3.3  (see  for  example  figure  13c).  When  simplified,  the  x  and 
y  terms  (first  and  second  terms  of  the  vector,  respectively)  remain  dependent  on  x. 
However,  the  z  term  (third  term  in  the  vector)  simplifies  to 

j,  =  n[ABc(xM  ~x0)  +  n%c(yM  ~  Vo)  +  n%c{zM  -  zp) 

\]  ( nABC )2  +  ( nABC )2  +  (nABc)2 


which  is  not  dependent  on  the  points  x  G  Sn-  As  a  result,  the  plane  containing 
the  points  x"  is  indeed  parallel  with  the  (x,y)  plane.  Further,  the  two-dimensional 
coordinate  system  depends  on  the  x  and  y  terms  of  x"  by 


x  = 


1  0  0 
0  1  0 

I*  ABC  •  (XM  —  Xp) 
n ABC  •  (x  -  Xo) 


1  0  0 
0  1  0 


M12(x  -  x0) 


where  x7  G  M2,  as  in  hgure  13d. 
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Recall  that  in  the  two-dimensional  coordinate  system,  the  RBF  interpolant  is 
defined  based  on  both  (f>(r],£,)  and  polynomial  terms  {717(77, in  (9).  Because 
some  x7  points  may  be  located  far  from  the  (0,  0)  (large  angle  between  thabc  and 
x  —  Xq )  the  polynomial  function  values  717(77,  £)  can  be  large  related  to  the  size  of  </>. 
As  a  result,  the  interpolation  matrix  A  for  s(r],£)  may  be  ill-conditioned.  To  remedy 
large  polynomial  terms,  the  points  x'  are  shifted  one  last  time  by  x'M,  the  projected 
midpoint  of  triangle  Tfe.  This  point  is  chosen  since  the  projected  neighboring  nodes 
are  nearest  Xm-  That  is,  the  final  two-dimensional  coordinate  system  is 


/  / 

=  x  ~XM 

A/r  ( nABC  ■  (XM  -  xo)  /  \  /  \ 

M12  - 7 - t-(x  -  Xo)  -  (XM  ~  x0) 

V  n ABC  ■  (x  -  Xo) 

1 

A/12 - 7 - C {jfi-ABC  X  ((x  -  Xo)  x  (xM  -  Xo))) 

TA-ABC  ■  (x  —  Xo) 

where  Xm  is  moved  to  the  origin.  Note  that  all  the  operations  performed  to  get  x  to 
lo  are  transformations  for  rotation  and/or  translation.  As  a  result,  the  planar  region 
Qfc  maintains  the  same  lengths  and  area  as  it  did  in  the  triangulation  T.  In  this 
process,  it  can  be  thought  of  that  the  surface  S  is  rotated  and  translated  (along  with 
the  integrand)  before  any  approximations  for  the  integral  (13)  are  carried  out.  With 
this,  the  planar  quadrature  weights  in  (13)  are  computed  identically  as  in 

section  2.2.3  with  calculations  of  RBF  integrals  evaluated  as  in  section  2.3.1. 


1  0  0 
0  1  0 
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(a) 


(d) 

Figure  13.  An  example  of  the  progression  of  coordinates  on  S  from  M3  in  Snk  to  M2  in  1\,. 

(a)  A  triangulation  of  the  surface 

h(x,  y ,  z)  =  ( x 2  +  y2  +  z2)2  —  2a2(x2  —  y2  —  z2)  +  a4  —  ft4  =  0  with  parameters  a  «  0.331, 
b  ~  0.348.  (b)  The  projection  of  J\fn  neighboring  nodes  and  the  corresponding  S^k  surface 
region  from  (a),  (c)  The  J\fn  projected  points  after  having  been  shifted  and  rotated,  (d) 
The  two-dimensional  coordinate  system  of  the  planar  region  fife  before  the  final  shift. 
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3.5  Converting  Planar  Weights  to  Surface  Weights 

Just  as  in  section  2.3.2,  this  section  will  take  the  surface  integral 

JJ  f(x,y,z)dS 
s*k 

for  Sak  the  region  of  interest  (i.e.  Snk  includes  the  neighboring  nodes  Afn  and  the 
corresponding  triangulation)  and  rewrite  it  as  an  area  integral  to  be  approximated 
with  the  quadrature  weights  determined  in  the  section  3.4.  The  conversion  from  a 
surface  to  area  integral  used  in  section  2.3.2  for  the  sphere  was  relatively  simple  as 
the  distance  between  the  projection  point  (origin)  and  the  projected  triangle  mid¬ 
point  xA,/  was  always  p,  and  the  nodes  on  S 'n  were  projected  onto  the  tangent  plane 
through  a  scalar  c  from  section  2.2.1.  For  an  arbitrary  smooth  closed  surface  S,  the 
derivation  for  the  transformation  ratio  for  p(p,£)  to  f(x,y,z )  focuses  on  the  use  of  a 
local  parameterization  of  S  described  by  x(p,£):  {x(r],l;),y(r),£),z(r},£)).  While  this 
chapter  acquires  much  of  its  information  from  [42],  this  particular  section  is  based  on 
the  details  in  [44], 

To  begin,  consider  the  points  x(p,  £),  x(p  +  A//,  £),  and  x(p,  £  +  A£)  on  S,  and  the 
corresponding  vectors  x(p+A rj,  £)—  x(p,  £)  and  x(p,  £+A £)— x(p,  £).  Let  A rj,  A^  >  0. 
Recall  that  when  two  vectors  a  and  b  compose  the  sides  of  a  parallelogram,  then  the 
area  of  the  parallelogram  is  A  =  ||a  x  b||.  Thus,  the  area  of  the  parallelogram  with 
sides  x(p  +  Ay,  f )  -  x(p,  0  and  x(p,  ^  +  A^)  -  x(p,  0  is 

As  patch  =  ||(x(r/  + Ap,0  -x(p,0)  X  (x(77,f  +  AO  -x(77,0)  l|2 
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which  can  be  rewritten  as 


A 


S  patch 


+  A'J]’  ^  -  x(^’  £))  x  W7?’  £  +  A0 

Ar)A£ 

x(r/  + A/7,0  -x(r/,Q  x  x(-/7,^  +  AQ  ~x(r/,Q 
Ar]  A£ 


-x(h,0)] 

A?/A£ 


2 


since  A 77,  A£  >  0.  As  Ary,  A£  — »  0, 


(h,0  x 


|x(,,0 


d//d£ 


where  dS1  is  the  same  used  for 


f(x,y,z)dS 


Sn,, 


In  planar  coordinates,  considering  the  same  parallelogram,  the  area  element  for  inte¬ 
gration  in  the  plane  is  given  as 


dA 


JU(r/,0  x  ^x(h,0 


drjd^ 


That  is, 


g(v,t,)dA 


Q,k 


such  that  <7(77,  £)  =  f(r},£)dS/dA  since  dS  =  ( dS/dA)dA .  Note  that  x  from  section 
3.3  is  used  instead  of  the  points  in  the  (77,  £)  plane.  As  described  before,  the  operations 
performed  in  section  3.4  preserve  the  same  area  between  planar  region  and  the 
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triangle  in  T.  Expanding  the  terms  in  dA, 


d  c\  3  UABC  ■  -  X0 

?)  =  Q-  (  Xo  + - — C\ - x  (XW  ?)  -  Xo) 

drj  dr)  \  nABc  ■  (x(ry,  0  -  x0) 


IlABC  •  (XM  ~  Xp) 

(nABC-  (x(r/,0  -x0))2 
-  ^nn_BC  •  /^(//-O)  (x(r/,0  -  xo) 


a 


(nABC  •  (x(ry,£)  -xo))^x(r/,0 


and  similarly  for  J|X(t7,£).  Therefore, 


|^x(7?,0  x  ^x(r/,0 


(nABC  •  (xm  ~  xq))2 
(nABC  •  (x(7/,0  -  x0))3 


(X(?Z,  0  -x0) 


^X(??,0  X  ^x(r/,0  )  )  nABC 


Substituting  this  back  into  the  equation  for  dA , 


JU(rz,0  x  ^x(r/,0 


dA  = 

1 1  (Jr)  '  ’ '  J '  d 

I nABC  •  (xm  —  Xq) | 2 

|nAsc-  (x(rz,0  -xo)|3 


dr]d£ 


(x(d,0  -  xo)  • 


|^x(rz,0  x  ^x(r/,0 


|nABc||2drzd^ 


It  was  noted  earlier  that  g(r),£)  =  f(r),£)dS/dA  since  dS  =  ( dS/dA)dA .  Hence, 


dS 


dS 

—dA 

dA 


£j*{v,0  x  ^x(?7,0 

drjd £ 

2 

lnABC"(xM— XC>)|2 

(xfa,0  -  Xo)  •  (^x(r?,0 

x  ix(^o) 

InAucIMr/d^ 

|nABC'(x(»?.t)-xo)|3 

nABC  /"y  y-T 

I 

(  n ABC  •  (x  -  Xo)  \ 

lln.4Bcl|2  ' 

§jx(vA)x-§tx(vA)  /  \ 

l^x(t/4)x^x(^)||2  v 

1 

\nABc  •  (xM  -  x0)/ 
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Therefore, 


xsnk  (/)  =  II  f(x,y,z)dS  = 

s^k 


1 use  /V 

-  Xo) 

1 f  f(x,y,z)  -j 

i|“ASc||2  ' 

/  nAsc  •  (x  -  x0) 

.8  x 
drj 

:(»7>0x^x(r?,0 

VnASC  '  (XM  -  X0) 

||^x(?7^)x^x(r?,0||2 


2 

dA  (29) 


Equation  (29)  is  then  used  for  a  parameterized  surface  S  defined  by  x(r/,  £).  From 
section  3.4, 


/(<?) 


g(Tl,  Qdridt 


J2wfBF^Vi.Q 

j= 1 


where  this  section  determined 


g(x,y,z)  =  f(x,y,z ) 


nASC  .  (V  y_'\ 

1|“ABC||S  ' 

|x(i/,^|x(,,0 

Xo) 

CM 

X 

uy> 

3 

*>1=0 

/  tlABC  •  (x  -  Xp) 
\nABC  •  (xM  —  Xo) 


As  a  result,  a  similar  procedure  as  in  section  2.3.2  applies  when  determining  the 
overall  weights  Wat.  Recall  that  each  quadrature  weight  wBBt  is  actually  computed 
as  ( wBBF)j  for  each  triangle  k  =  1,. and  each  projected  node  ( Xj,yj,Zj )  for 
j  —  As  with  the  spherical  quadrature,  let  /C?; ,  for  i  =  quadrature 

nodes,  be  the  set  of  all  pairs  (k,j)  such  that  {(Vk)j,  (£k)j)  l— ^  ( x^y^Zi ).  Then  the 


57 


integral  of  /  over  the  whole  surface  S  is 


w)=E// 

k= i 


f(x,  y,  z)dS 


E 


E 


v 


(k,j)eiCi 


.  /V-y-,  ^  Y  j  ^ 

IlnAsdb  dX*b  lXOjfcj 

£jx(v,Ox-§e*(ihO  ,Ur  x 

:o)fc) 

f  n  \bc  ■  ((xfc)j  ^  (x0)fc) 
\nABC  •  ((xjf)t  —  (xo)fc) 


A 

/ 


f(xi,yi,Zi ) 

N 

=  Y  2/*>  Zi )  =  ^S(/) 
2=1 


(30) 


Equation  (30)  is  the  numerical  approximation  for  the  integral  of  /  over  an  explicitly 
parameterized  surface.  Suppose  instead  that  the  surface  S  is  defined  implicitly  by 
h(x,y,z )  =  0.  Note  that 

£j*(v,0  x  ^xfaO 

ll|^x(h,0  X  J|X(77,0 lb 

is  a  unit  normal  to  the  surface  S'.  This  can  be  related  to  the  gradient  S7h(x,y,  z) 
(also  a  surface  normal)  by 


i^xObO  x  ^x(h,Oll2 

=  sign 


(h,0  x 


S7h(x,  y,  z] 


Vh(x,y,z) 

\\Vh(x,y,z)\\2 


which  when  placed  in  (30)  gives 


Wf)  =  tjj 

k= 1 '' 


f{x,  y,  z)dS 


N 


E 

i=  1 


E  «BF), 

(k,j)eiCi 


iS  •  ~ 

v/f(w)||2  •  ((x^)i  -  (xo)fc) 


f  n ABC  ■  ((Xfc)j  -  (xp)fc) 
\nABC  ’  ((xM)fc  ~  (Xo)fe) 


•  f(xi,yi,Zi ) 

N 

=  X  Vi’  Zi )  =  ^s(/) 

Z=1 


(31) 


Equation  (31)  is  hence  the  numerical  approximation  for  the  integral  of  /  over  an 
implicitly  defined  surface.  Consider,  however,  that  S  may  not  be  defined  explicitly 
through  parameterization  nor  implicitly  by  h(x,y,z )  =  0.  Just  as  the  intermediate 
quadrature  weights  in  section  3.3  were  determined  with  only  quadrature  nodes  Sn 
and  a  triangulation  T,  Is(f)  must  also  be  computable  under  the  same  conditions. 
Section  3.6  covers  this  case. 


3.6  Approximating  a  Normal  to  the  Surface 

Suppose  that  the  surface  S  is  defined  neither  explicitly  through  a  parameterization 
nor  implicitly  by  h(x,  y ,  z)  =  0.  Instead,  only  a  node  set  Sn  and  triangulation  T  are 
provided.  The  terms 

Vh(x)  |%Ox|x('f;,0 

Si an  ni^x(//,o  x  ^x(//,oii2 

were  the  only  terms  in  (30)  and  (31)  which  did  not  use  the  information  provided  by 
Sn  and  T.  Notice  that  just  as  points  in  the  (?/,£)  plane  are  dependent  on  the  nodes 
in  Sn,  so  too  can  the  nodes  {(xj,yj,  Zj)}™=1  for  n  projected  nodes  be  expressed  as 
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a  function  of  77  and  £.  That  is,  the  projection  and  transformations  carried  out  on 
{( Xj ,  y3 ,  Zj)}"=1  provide  a  local  “explicit”  parameterization  for  each  projected  point  in 

Sn- 

ix(Vj,Q,y(Vj,Q,z(Vj,Q)  =  for  j  =  1 

where  (r)j,£j)  is  the  two-dimensional  planar  version  of  ( Xj,yj,Zj ).  While  the  true 
parameterization  (x(r),  £),  y(rj,  £),  z(rj,  £))  is  unknown  (since  S  is  not  equation-based), 
(x(rj,  £),  y(r),  £),  z(r),  £))  can  be  approximated  through  a  radial  basis  function  interpo¬ 
lation 


n  /  \  M 

x(v,0  ~  sx(r],0  :=  J2c%3F<f>  (  ^ (v  -  Vj)2  +  (£  -  £j)2) 

j= 1  '  '  1=1 

n  /  \  M 

y(v,0  ■=  Y  cyfF $  ^')2  +  (f  ~  0)2 )  +'%2<%Mv,0 

j= 1  ^  ’  1=1 

n  /  \  M 

z(v, 0  ~  sz{n, 0  -.=  Y czfF$  ( \J (v -  vj)2 ~+if- 0)2 )  +  Y ^zMv, 0 

j= i  ^  '  1=1 

Without  loss  of  generality,  cFfF,...,cFFF,  c^1,...,c^M  G  M  are  chosen  such  that 
sx(r}j,€j)  =  x(Vj^j),  3  =  1,2 along  with  constraints  V;  icxfF'Mrlj,€j)  =  0, 
for  l  —  1,2,  The  coefficients  for  sx(r),£),  sy(r},^),  and  sz(r},£)  are  evaluated  by 

the  matrix  multiplication  elaborated  in  section  2.2.3.  Once  the  interpolants  sx(rj,^), 
sy(j 7,0,  an<l  sz(y,0  are  constructed, 


9 


X 


can  be  approximated  by 


d_ 

drj 


s*{v,£)  x 


d 
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where 


Sx  (v,0 

x(v,0 

sx(v,Q  = 

Sy(v,0 

r^i 

y(v,0 

Sz(V,0 

z(v,0 

The  partial  derivatives  of  sx(r),£,)  (and  therefore  <sx (?/,£))  can  be  determined  by,  for 
instance, 


r\  Tl  r\  /  , - \  M  r\ 

=  E>jfF— </.  (y  (v  -  vi)2 + (e  -  e,)2j  + 

q^s-(V,0  =  ^4^—0  (y (b  -  hi)2  +  (f  “  0)2J  +  EC^7r^’£) 

where  the  partial  derivatives  of  the  RBF  0  will  depend  on  the  chosen  basis  function. 
Hence,  (30)  becomes 


A' 


W)  =  E 


k= 1 


f(x,  y,  z)dS 


E 


/ 

E 


Wn^ch  ■  ^ 

1 

(  n ABC  ■  ((Xfc)j  -  (xo)fc) 

\\&jS*(y,t)x%sAv,t) lb  1 o)k) 

1 

\nABC  '  ((xM)fc  —  (x0)fc) 

\ 

/ 


f(xi,  Vi,  Zi) 

N 

—  E  Wif(Xi,  Vi,  =  is(f) 

i= 1 


(32) 


which  is  the  numerical  approximation  of  /  over  a  surface  S'  expressed  by  a  set  of 
nodes  S>at  and  triangulation  T ■ 
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IV.  Results 


4.1  Test  Surface 

This  chapter  applies  the  material  presented  in  Chapter  III  to  a  surface  S  implicitly 
defined  by 


h(x,  y,  z )  =  ( x 2  +  y2  +  z2)2  —  2 a2(x2  —  y2  —  z2)  +  a4  —  b 4  =  0 


(33) 


and  several  test  integrands.  S'  is  a  surface  of  revolution  induced  by  the  Cassini  oval 
[45]  in  the  (x,  y)  plane  (rotated  about  the  rc-axis)  with  parameters  a  and  b.  Examples 
can  be  seen  in  figure  15.  Parameters  a  and  b  are  chosen  to  satisfy  a  =  X b  for  values 
of  A  =  (0.8,  0.85,  0.9,  0.95}  and  b  such  that  S  has  surface  area  equal  to  1  (see  figure 
14  for  slices  of  S  in  the  (x,y)  plane).  This  was  done  to  normalize  the  surface  area’s 
contribution  to  any  error  curves  shown.  The  surface  area  of  S  can  be  approximated 
using  MATLAB’s  quad2d  where,  in  polar  coordinates  [39], 


Ix(W)x|x(W) 


d0d(f> 

2 


with 


and  [45] 


x(M) 


r(<t>) 


COS  (4>) 

sin(0)sin(6)) 
sin(4>)cos(9 ) 


r(0)  =  yj  \[¥ 


a4  +  a4cos2(20)  +  a2  cos(20) 


The  quadrature  nodes  Sn  and  triangulation  T  were  generated  using  distmeshsur- 
face  [43],  where  the  Ell  distance  (as  described  in  section  2.3.2)  varied  from  0.005  to 
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X 

Figure  14.  Slices  of  Cassini  ovals  used  for  application,  where  a  =  Xb  for  A  values  shown. 
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a  =  0.95  and  6  =  1 


a  =  0.6  and  6  =  1 


a  =  1.05  and  6  =  1 


Figure  15.  Example  Cassini  ovals  rotated  about  the  x-axis  for  given  parameter  values  a 
and  b.  For  a  <  b,  S  appears  as  an  increasingly  pinched  oval  as  a  — >  b.  For  a  =  b,  S  is 
completely  pinched  (but  inherently  not  smooth),  while  a  >  b  depicts  S  as  a  piecewise 
smooth  closed  surface  with  two  separated  ovals. 

0.10.  Mesh  sizes  much  smaller  than  0.005  were  not  considered  because  of  the  com¬ 
putational  time  required  to  generate  them  with  distmeshsurface,  while  mesh  sizes 
much  larger  than  0.10  generate  too  few  nodes  for  a  reasonable  approximation. 

4.2  Applying  the  Chapter  III  Quadrature  Method 

Following  section  3.4,  the  approximation  (12)  here  uses  the  RBF  </>(r)  =  r7,  n  =  80 
(number  of  nearest  neighbors),  and  m  —  7  (degree  of  polynomial  terms).  From  the 
implicitly  defined  surface  in  (33),  the  quadrature  weights  for  (31)  were  determined. 
If  f(x,y,z)  =  1,  then  (31)  becomes  the  surface  area  of  S\ 


As  described  in  section  4.1,  the  surface  area  of  S  was  set  up  to  equal  1.  With  this,  the 
sum  of  the  quadrature  weights  in  (33)  should  be  approximately  equal  to  1.  Suppose, 
however,  that  (32)  was  used  instead  such  that  an  approximation  to  the  true  normal 
(' Vh(x,y,z ))  was  utilized: 
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Error  in  Quadrature  Over  S  with 
Exact  Surface  Normal:  Testing  Surface  Area 


logio(lV) 

(a) 


Error  in  Quadrature  Over  S  with 
Approximated  Surface  Normal:  Testing  Surface  Area 


logic  (AQ 
(b) 

Figure  16.  Relative  error  in  RBF  quadrature  for  the  surface  area  over  surfaces 
h(x,  y,  z)  =  0  for  various  A  on  sets  of  quasi-uniformly  spaced  nodes  with:  (a)  surface 
normal  computed  via  V/i;  (b)  surface  normal  approximated.  For  the  computations, 

(j){r)  =  r7,  n  =  80,  and  m  =  7. 
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Estimates  for  computing  surface  areas  of  S  defined  by  (33)  are  given  in  figure  16.  In 
either  case,  a  convergence  rate  of  0(N~ 3'5)  is  achieved.  The  large  errors  for  smaller 
N  in  both  cases,  and  larger  N  in  the  case  of  the  approximate  surface  normal,  will  be 
discussed  after  the  results  are  presented. 

Now  consider 

fi(x,y,z)  =  F,-  •  n,i  =  1,2,3 

where 


n  = 


Vh(x,y,z) 

\\Vh(x,y,z)\\2 

x 


Fi  =  - 

3 


Fo  = 


(. x 2  +  y2  +  z2y/2 


F 3  =  V  x  F36  =  V  x 


-\y 


\x 


Note  that  f\ ,  /3  G  C°°(S )  while  /2  is  also  continuous  except  for  when  the  denominator 
is  zero.  That  is,  /2  G  C°°(S)  since  x2  +  y2  +  z2  =  0  does  not  exist  on  S.  Also  consider 


h{x,y,z) 


2 

— tan_1(100^)  and  f5(x,y,z)  =  sign(^)  =  { 
7 r 


+1; 

0, 

-1, 


z  >  0 
z  =  0 
z  <0 
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Mathematically,  X s{fi)  is  (by  the  Divergence  Theorem)  [39] 


ls(fi)=  J(F1.n)dS=  IIJ  (V-F ,)dV 
s  v 

•  r 

1  dV 

! 

v 

That  is,  it  is  the  volume  of  the  region  enclosed  by  the  surface  S.  Recall  that  S  is  a 
surface  of  revolution  about  the  rr-axis,  such  that 


h(x,  y,  0)  —  (x  +y  —  2a  (x  —  y  )  +  a  —6=0 


(34) 


Then, 


r{x )  =  y(x )  =  \j —a?  —  x2  +  VbA  —  4  a2x2 


where  y{x)  is  one  of  four  roots  of  (34)  and,  from  this,  x  =  \/  a2  +  62  when  y  =  0. 
Therefore, 


^(/i)  = 


1 dydz  dx 


x  \  A 
V  a2+62 


7r r2{x)dx 


\/  a2+62 


=  2vr 

7T 

6  a 


a2  —  x2  +  -\/64  —  4a2x2j  dx 


2a (62  —  2a2)  V a2  +  62  +  364sinh  1  ( 

0“ 


The  error  plot  for  Xg(/i)  can  be  found  in  figure  17. 

The  integral  of  /2(:r,  y,  z ),  however,  uses  the  Divergence  Theorem  [39]  and  a  prop- 
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Error  in  Quadrature  Over  S  with 
Exact  Surface  Normal:  Test  Function  /i 


Error  in  Quadrature  Over  S  with 
Approximate  Surface  Normal:  Test  Function  f\ 


logic  (N) 

(b) 

Figure  17.  Relative  error  in  RBF  quadrature  for  function  f\(x,y,z)  over  surfaces 
h(x,  y,z)  =  0  for  various  A  on  sets  of  quasi-uniformly  spaced  nodes  with:  (a)  surface 
normal  computed  via  V/i;  (b)  surface  normal  approximated.  For  the  computations, 

(j>{r)  =  r7,  n  =  80,  and  m  =  7. 


erty  concerning  F 2(x,y,z)  [46]: 


Xs(/2)  =  If  (F2  •  n)  dS  =  HI  (V  ■  F2)  dV 
s  v 

=  HI  And3  ( x 2  +  y2  +  ^2)  dV 
v 

=  An 


The  error  plot  for  I$'(/2)  can  be  found  in  figure  18. 

The  integral  of  fs(x,y,z)  uses  Stokes’  Theorem  [39], 

Mh)  =  If  ((V  X  F3b)  •  n)  dS  =  I F36  •  dr 

S  dS 

tf 

=  I  F36(r (t))  ■  r'(t)dt 

to 

tf 

to 

=  I  (f3^}  [x,  y,  z)dx  +  F3(b2)  [x,  y,  z)dy ) 
as 

which  Green’s  Theorem  [39]  relates  the  last  line  to 

xsih)  =  If  (J^Fg\x,y,z)  -  y,  z)^j  dxdy 

A 


dt 


The  last  step  follows  from  the  assumptions  of  Stokes’  Theorem,  where  the  the  surface 
contour  (boundary)  traverses  one  direction  for  the  upper  half  of  S  but  the  opposite 
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Error  in  Quadrature  Over  S  with 
Exact  Surface  Normal:  Test  Function  /2 


Error  in  Quadrature  Over  S  with 
Approximate  Surface  Normal:  Test  Function  f  > 


logic  (N) 

(b) 

Figure  18.  Relative  error  in  RBF  quadrature  for  function  f2(x,y,z)  over  surfaces 
h(x,  y,z)  =  0  for  various  A  on  sets  of  quasi-uniformly  spaced  nodes  with:  (a)  surface 
normal  computed  via  V/i;  (b)  surface  normal  approximated.  For  the  computations, 

(j>{r)  =  r7,  n  =  80,  and  m  =  7. 
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Error  in  Quadrature  Over  S  with 
Exact  Surface  Normal:  Test  Function  /3 


Error  in  Quadrature  Over  S  with 
Approximate  Surface  Normal:  Test  Function  /3 


logio(lV) 

(b) 

Figure  19.  Relative  error  in  RBF  quadrature  for  function  fz(x,y,z)  over  surfaces 
h(x,  y,z)  =  0  for  various  A  on  sets  of  quasi-uniformly  spaced  nodes  with:  (a)  surface 
normal  computed  via  V/i;  (b)  surface  normal  approximated.  For  the  computations, 

4>(r)  =  r7,  n  =  80,  and  m  =  7. 
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direction  for  the  bottom  half  of  S',  so  that  the  surface  normals  are  consistent  from 
one  to  the  other  [39] .  This  causes  the  integrals  for  the  upper  and  lower  halves  of  S  to 
be  equal  in  magnitude  but  opposite  in  sign.  The  resultant  error  plot  for  can 

be  found  in  figure  19. 

Now,  consider  f4(x,y,z)  =  ^tan_1(100^)  (where  the  details  concerning  are 
from  [47]).  Since  S'  is  a  surface  of  revolution  induced  by  the  Cassini  oval  in  the  (x,y) 
plane,  then  S'  G  M3  is  also  symmetric  about  the  (x,  y )  plane.  Therefore,  define  the 
portion  of  S  for  which  z  >  0  to  be  S+ .  Similarly,  let  S'-  be  the  portion  of  S  for  which 
z  <  0.  So, 


MU)  = 


—tan  1(100^)rfS 

7T 


S+ 


—tan  1(100z)<iS' +  J J  —tan  1(100z)dS' 
s- 


— tan  1(100z)(iS' +  JJ  —  tan  x(100( — ^))c?Sf 
s+  s+ 


Since  arctangent  is  an  odd  function  (i.e.  tan  :(— x)  =  —tan  1(x)), 


ls(f4)  =  JJ  —tan  1(100z)dS  —  JJ  —tan  1(100x;)(iS  =  0 

s+  s+ 


The  resultant  error  plot  for  f4(x,y,z )  is  illustrated  in  figure  20. 
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Error  in  Quadrature  Over  S  with 
Exact  Surface  Normal:  Test  Function  /4 


Error  in  Quadrature  Over  S  with 
Approximated  Surface  Normal:  Test  Function  /4 


logio(lV) 

(b) 

Figure  20.  Relative  error  in  RBF  quadrature  for  function  f4(x,y,z)  over  surfaces 
h(x,  y,z)  =  0  for  various  A  on  sets  of  quasi-uniformly  spaced  nodes  with:  (a)  surface 
normal  computed  via  V/i;  (b)  surface  normal  approximated.  For  the  computations, 

4>(r)  =  r7,  n  =  80,  and  m  =  7. 
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Similar  to  X5(/4),  the  integral  of  f5(x,y,z)  =  sign(z)  is 


Zs{h)  =  jj  sign (z)d,S 
s 

=  IJ  (+l)dS  +  JJ(- l)dS 

s+  s- 


s+  s- 


Since  the  surface  integral  over  S  (with  an  integrand  of  1)  represents  the  surface  area 
of  S,  the  surface  integrals  over  S+  and  S~  each  represent  half  the  surface  area  of  S. 
That  is, 


The  resultant  error  plot  for  f${x,y,z)  is  illustrated  in  figure  21. 

Notice  that  with  almost  all  the  error  plots,  the  convergence  rate  is  0(N~ 3'5)  or  bet¬ 
ter,  which  is  that  attained  with  the  sphere  for  continuous  functions.  The  convergence 
rate  for  f±(x,y,z)  is  closer  to  0(N~2-5).  A  lower  convergence  rate  is  to  be  expected 
due  to  the  steep  gradient  that  occurs  along  the  (x,  y)  plane.  The  convergence  rate 
for  fs(x,y,z)  was  approximately  O(N~075),  which  is  similar  to  the  convergence  rate 
illustrated  by  f2(x,y,z)  in  section  2.3.3.  These  rates  are  not  surprising,  since  as  long 
as  the  surface  gradient  V h  is  nonzero,  the  surface  S  is  as  smooth  as  the  sphere. 

Recall  the  computational  costs  of  O(NlogN)  operations  and  O(N)  memory  stor¬ 
age  for  the  sphere  quadrature.  These  computational  costs  will  remain  nearly  the  same 
for  the  smooth  closed  surface  S  since  the  operations  performed  for  the  intermediate 
quadrature  weights  in  (13)  and  final  quadrature  weights  in  (1)  are  the  same  for  both 
the  sphere  and  smooth  surface  S.  The  cases  of  the  smooth  surface  and  the  sphere 
differ  only  in  that  a  projection  origin  must  be  computed  for  each  triangle.  This  adds 
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Error  in  Quadrature  Over  S  with 
Exact  Surface  Normal:  Test  Function  /5 


logio(JV) 

(a) 


Error  in  Quadrature  Over  S  with 
Approximated  Surface  Normal:  Test  Function  f5 


logic  (AQ 
(b) 

Figure  21.  Relative  error  in  RBF  quadrature  for  function  f$(x,y,z)  over  surfaces 
h(x,  y,z)  =  0  for  various  A  on  sets  of  quasi-uniformly  spaced  nodes  with:  (a)  surface 
normal  computed  via  V/i;  (b)  surface  normal  approximated.  For  the  computations, 

4>(r)  =  r7,  n  =  80,  and  m  =  7. 
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only  a  few  scalar  operations  and  an  initial  O(NlogN)  sorting  process  (to  determine 
the  vectors  n ab,  n bc,  and  ncA  hi  section  3.2). 

For  fixed  m  —  7,  n  —  80,  and  0(r)  =  r  ,  the  total  cost  in  computing  the  RBF 
quadrature  weights  over  S  implicitly  defined  by  (33)  is  illustrated  in  figure  22.  Since 
the  spherical  quadrature  method  and  the  Chapter  111  quadrature  method  differ  only 
by  the  computation  of  the  projection  origin,  it  is  not  surprising  that  they  exhibit  sim¬ 
ilar  computational  costs.  Whether  using  the  exact  gradient  or  approximate  to  it,  the 
Chapter  111  quadrature  method  is  computationally  inexpensive  at  O(N)  operations 
and  0(N )  memory  usage.  The  memory  usage  plot  exhibits  no  apparent  difference 
between  the  two,  while  the  computation  time  plot  exhibits  a  constant  displacement. 
This  displacement  comes  from  computing  the  gradient  approximation.  While  parallel 
scaling  was  not  tested  for  the  method  in  Chapter  111,  it  is  expected  parallel  scalabil¬ 
ity  with  number  of  cores  will  be  observed  similar  to  the  spherical  quadrature  plot  in 
figure  8c. 

In  each  of  the  error  plots,  the  error  is  large  for  N  less  than  roughly  102'6.  This  is 
because  mesh  sizes  that  are  too  large  generate  too  few  points  (N  ~  102  6  ~  400)  such 
that  the  requirement  for  Gnomonomic  projection  (x-xo  within  90°  of  n abc)  is  often 
violated  for  n  =  80.  Meanwhile,  increases  in  error  past  N  «  104'5  ~  31,  500  (in  the 
case  of  the  approximate  surface  normal)  are  indicative  of  numerical  approximation 
errors  when  computing  derivatives  for  (32).  The  deliberation  for  the  increase  in  error 
is  left  for  future  consideration. 
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Figure  22.  Timing  results  for  computation  of  the  RBF  quadrature  weights  for  the  surface 
area  of  S  implicitly  defined  by  (33)  as  applied  in  section  4.2.  The  quadrature  method  has 
rate  O(N)  for  computation  cost  and  memory  usage  for  at  least  millions  of  nodes.  An 
additional  computational  cost  is  observed  when  computing  the  approximation  to  the 
gradient.  These  plots  were  generated  in  machines  with  dual  Intel  Xeon  E5-2687W  3.1 

GHz,  8-core  processors. 
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V.  Conclusion 


5.1  Future  Considerations 

With  such  rich  areas  of  applicability,  this  thesis  infers  many  avenues  for  future 
consideration.  The  work  here  uses  only  monomial  radial  basis  functions,  whereas 
infinitely  smooth  RBFs  have  more  robust  convergence  properties  [7].  Therefore,  it 
is  important  to  consider  the  stable  use  of  these  RBFs  in  the  future,  as  was  done  for 
interpolation  in  [13,  14,  48,  49,  50]. These  infinitely  smooth  RBFs  include  the  shape 
parameter  e  that  has  been  discussed  throughout  the  literature  [49,  51,  52],  Using 
different  RBFs,  and  analyzing  how  £  affects  the  results  would  provide  intriguing  future 
work. 

Decreasing  computation  time  when  generating  either  Sn  or  T  is  a  separate  but 
important  issue.  Using  a  less  costly  version  of  distmeshsurface  would  help  with 
computing  these.  Furthermore,  distmeshsurface  does  not  generate  Sn  and  T  for  a 
surface  defined  by  parameterized  coordinates  x(u,  v ),  y(u,  v),  and  z(u,  v).  Developing 
an  algorithm  that  does  so  could  be  a  significant  contribution. 

When  analyzing  the  numerical  approximation  for  errors  for  an  approximated  sur¬ 
face  normal  (see  section  4.2),  there  is  an  increase  in  error  past  N  «  104'5  ~  31,500. 
Deliberating  the  source  of  this  is  an  important  avenue  for  future  research. 

Other  future  considerations  include  non-smooth  (jagged)  or  non-closed  surfaces. 
Non-smooth  surfaces  pose  an  issue  when  approximating  the  derivative  in  use  with 
(32).  Furthermore,  the  quadrature  method  from  this  thesis  requires  the  projection 
of  nearest  neighbors,  whereby  cusps  in  the  surface  may  cause  difficulty.  Non-closed 
surfaces  present  a  similar  issue  when  using  points  near  the  surface  boundary.  Approx¬ 
imating  the  integral  of  a  scalar  function  /  over  a  non-smooth  or  non-closed  surface 
would  indeed  be  a  challenging  but  significant  endeavor. 
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5.2  Concluding  Remarks 


In  the  case  of  approximating  surface  integrals  on  the  sphere,  the  convergence  rate 
of  the  method  in  Chapter  II  was  determined  to  be  on  the  order  of  0(N~3'5)  for  quasi¬ 
uniform  node  sets.  Therefore,  the  spherical  quadrature  method  from  Chapter  II  yields 
accurate  results  with  fast  convergence.  The  results  of  Chapter  IV  showed  simlar 
convergence  for  smooth,  closed  surfaces  at  still  low  computational  cost.  However, 
the  computation  for  generating  the  quadrature  nodes  Sn  and  triangulation  T  via 
distmeshsurface  took  an  intensive  amount  of  time.  It  is  only  after  these  nodes  and 
triangulation  are  generated  that  the  results  in  Chapter  IV  can  be  interpreted. 

The  method  for  using  an  approximation  to  the  surface  normal  is  very  robust  in 
the  sense  that  the  integral  of  a  scalar  function  /  should  be  able  to  be  approximated 
over  any  smooth  closed  surface  S  regardless  of  a  (un)known  parameterization.  As 
the  accuracy  of  the  approximation  attains  levels  of  approximately  10-10,  this  can  be 
considered  an  excellent  approximation.  Furthermore,  the  convergence  rate  remained 
the  same  at  0(N~3'5).  As  the  method  discussed  in  Chapter  III  allows  the  surface 
to  be  on  scattered  grids,  it  can  be  used  for  a  wide  range  of  applications.  Depending 
on  the  surface,  however,  the  user  may  need  to  acquire  a  large  number  of  samples  to 
achieve  a  desired  accuracy. 


79 


Bibliography 


1.  R.  L.  Hardy,  “Theory  and  Applications  of  the  Multiquadric-Biharmonic  Method,” 
Comput.  Math.  Appl .,  vol.  19,  no.  8-9,  pp.  163-208,  1990. 

2.  J.  A.  Reeger  and  B.  Fornberg,  “Numerical  Quadrature  over  the  Surface  of  a 
Sphere,”  Stud.  Appl.  Math.,  2016.  DOI:  10. 1111/sapm.  12106. 

3.  Z.  P.  Bazant  and  B.  H.  Oh,  “Efficient  Numerical  Integration  on  the  Surface  of  a 
Sphere,”  Zamm-Z.  Angew.  Math.  ME.,  vol.  66,  pp.  37-49,  1986. 

4.  A.  Sommariva  and  R.  Womersley,  “Integration  by  RBF  over  the  Sphere,”  Appl. 
Math.  Report  AMR05/17 ,  2005. 

5.  R.  S.  Womersley  and  I.  H.  Sloan,  “Interpolation  and  Cubature  on  the  Sphere.” 
http://web.maths.unsw.edu.au/  rsw/Sphere/,  2003/2007. 

6.  J.  P.  Bardhan,  Efficient  Numerical  Algorithms  for  Surface  Formulations  of  Math¬ 
ematical  Models  for  Biomolecule  Analysis  and  Design.  PhD  thesis,  Massachusetts 
Institute  of  Technology,  2006. 

7.  M.  D.  Buhmann,  “Radial  Basis  Functions,”  Acta  Numerica,  vol.  9,  pp.  1-38,  2000. 

8.  J.  C.  Mairhuber,  “On  Haar’s  Theorem  Concerning  Chebychev  Approximation 
Problems  Having  Unique  Solutions,”  P.  Am.  Math.  Soc.,  vol.  7,  no.  4,  pp.  609- 
615,  1956. 

9.  B.  Fornberg,  E.  Larsson,  and  N.  Flyer,  “Stable  Computations  with  Gaussian 
Radial  Basis  Functions,”  SIAM  J.  Sci.  Comput.,  vol.  33,  pp.  869-892,  2011. 

10.  B.  Fornberg  and  J.  Zuev,  “The  Runge  Phenomenon  and  Spatially  Variable  Shape 
Parameters  in  RBF  Interpolation,”  Comput.  Math.  Appl,  vol.  54,  pp.  379-398, 
2007. 

11.  S.  Bochner,  “Monotone  Funktionen,  Stieltjessche  Integrate  und  Harmonische 
Analyse,”  Math.  Ann.,  vol.  108,  pp.  378-410,  1933. 

12.  I.  J.  Schoenberg,  “Metric  Spaces  and  Positive  Definite  Functions,”  T.  Am.  Math. 
Soc.,  vol.  44,  pp.  522-536,  Nov.  1938. 

13.  R.  L.  Hardy,  “Multiquadric  Equations  of  Topography  and  Other  Irregular  Sur¬ 
faces,”  J.  Geophys.  Res.,  vol.  76,  pp.  1905-1915,  March  1971. 

14.  E.  J.  Kansa,  “Multiquadrics  -  A  Scattered  Data  Approximation  Scheme  with 
Applications  to  Computational  Fluid-Dynamics  -  I:  Surface  Appoximations  and 
Partial  Derivative  Estimates,”  Comput.  Math.  Appl,  vol.  19,  no.  8-9,  pp.  127-145, 
1990. 


80 


15.  A.  I.  Tolstykh,  “On  Using  RBF-Based  Differencing  Formulas  for  Unstructured 
and  Mixed  Structured-Unstructured  Grid  Calculations,”  in  Proceedings  of  the 
16th  IMACS  World  Congress  on  Scientific  Computation,  Applied  Mathematics 
and  Simulation ,  p.  6,  2000.  Lausanne,  Switzerland. 

16.  J.  G.  Wang  and  G.  R.  Liu,  “A  Point  Interpolation  Meshless  Method  Based  on 
Radial  Basis  Functions,”  Int.  J.  Numer.  Meth.  Eng.,  vol.  54,  no.  11,  pp.  1623- 
1648,  2002. 

17.  G.  B.  Wright,  Radial  Basis  Function  Interpolation:  Numerical  and  Analytical 
Developments.  PhD  thesis,  LIniversity  of  Colorado  at  Boulder,  2003. 

18.  C.  Shu,  H.  Ding,  and  K.  S.  Yeo,  “Local  Radial  Basis  Function-Based  Differential 
Quadrature  Method  and  its  Application  to  Solve  Two-Dimensional  Incompressible 
Navier-Stokes  Equations,”  Comput.  Method  Appl.  M.,  vol.  192,  no.  7,  pp.  941- 
954,  2003. 

19.  N.  Flyer,  G.  B.  Wright,  and  B.  Fornberg,  “Radial  Basis  Function-Generated  Finite 
Differences:  A  Mesh-Free  Method  for  Computational  Geosciences,”  Handbook  of 
Geomathematics,  2014. 

20.  P.  Keast  and  J.  Diaz,  “Fully  Symmetric  Integration  Formulas  for  the  Surface  of 
the  Sphere  in  S  Dimensions,”  SIAM  J.  Numer.  Anal.,  vol.  20,  no.  2,  pp.  406-419, 
1983. 

21.  V.  Lebedev,  “Quadratures  on  a  Sphere,”  USSR  Comp.  Math.  Math.+,  vol.  16, 
pp.  1-23,  1976. 

22.  A.  Stroud,  Approximate  Calculation  of  Multiple  Integrals.  Prentice-Hall,  1971. 
New  Jersey,  USA. 

23.  B.  Fornberg  and  J.  Martel,  “On  Spherical  Harmonics  Based  Numerical  Quadra¬ 
ture  over  the  Surface  of  a  Sphere,”  Adv.  Comput.  Math.,  vol.  40,  no.  5-6,  pp.  1169- 
1184,  2014. 

24.  K.  E.  Atkinson,  “Numerical  Integration  on  the  Sphere,”  J.  Aust.  Math.  Soc.  B, 
vol.  23,  pp.  332-347,  1982. 

25.  E.  Fuselier,  T.  Hangelbroek,  F.  J.  Narcowich,  J.  D.  Ward,  and  G.  B.  Wright, 
“Kernel  Based  Quadrature  on  Spheres  and  Other  Homogeneous  Spaces,”  Numer. 
Math.,  vol.  127,  pp.  57-92,  2014. 

26.  D.  Chien,  “Numerical  Evaluation  of  Surface  Integrals  in  Three  Dimensions,” 
Math.  Comput.,  vol.  64,  no.  210,  pp.  727-743,  1995. 

27.  A.  Sommariva  and  M.  Vianello,  “Meshless  Cubature  by  Green’s  Formula,”  Appl. 
Math.  Comput .,  vol.  183,  pp.  1098-1107,  2006. 


81 


28.  J.  D’Elia,  L.  Battaglia,  A.  Cardona,  and  M.  Storti,  “Full  Numerical  Quadra¬ 
ture  of  Weakly  Singular  Double  Surface  Integrals  in  Galerkin  Boundary  Element 
Methods,”  Int.  J.  Numer.  Meth.  Eng.,  vol.  27,  no.  2,  pp.  314-334,  2011. 

29.  G.  Green,  An  Essay  on  the  Application  of  Mathematical  Analysis  to  the  Theories 
of  Electricity  and  Magnetism.  Author,  1828.  Nottingham. 

30.  A.  Gonzalez,  “Measurement  of  Areas  on  a  Sphere  using  Fibonacci  and  Latitude- 
Longitude  Lattices,”  Math.  Geosci .,  vol.  42,  no.  1,  pp.  49-64,  2010. 

31.  D.  P.  Hardin  and  E.  B.  Saff,  “Discretizing  Manifolds  via  Minimum  Energy  Points,” 
Notices  of  the  AMS,  vol.  51,  no.  10,  pp.  1186-1194,  2004. 

32.  S.  V.  Borodachov,  D.  P.  Hardin,  and  E.  B.  Saff,  “Low  Complexity  Methods  for 
Discretizing  Manifolds  via  Riesz  Energy  Minimization,”  Found.  Com.put.  Math., 
vol.  14,  no.  6,  pp.  1173-1208,  2014. 

33.  R.  J.  Renka,  “Algorithm  772:  STRIPACK:  Delaunay  Triangulation  and  Voronoi 
Diagram  on  the  Surface  of  a  Sphere,”  ACM  T.  Math.  Software,  vol.  23,  pp.  416- 
434,  Sept.  1997. 

34.  B.  Delaunay,  “Sur  la  Sphere  Vide,”  B.  Acad.  Sci.  USSR,  vol.  7,  pp.  793-800,  1934. 

35.  T.  M.  Liebling  and  L.  Pournin,  “Voronoi  Diagrams  and  Delaunay  Triangulations: 
Ubiquitous  Siamese  Twins,”  Doc.  Math.  Extra  Volume:  Optimization  Stores, 
pp.  419-431,  2012. 

36.  J.  P.  Snyder,  Map  Projections-A  Working  Manual,  vol.  1395.  US  Government 
Printing  Office,  1987.  Washington,  D.C.,  USA. 

37.  J.  A.  Reeger,  “Handout:  19  December  2013.”  Private  communication,  Dec.  2013. 

38.  J.  A.  Reeger  and  B.  Fornberg,  “Numerical  Quadrature  Over  the  Surface  of  a 
Sphere.”  Private  communication,  Nov.  2014. 

39.  J.  Stewart,  Calculus:  Early  Trans cendentals.  Brooks/Cole-Thomson  Learning, 
5  ed.,  2003.  Belmont,  CA,  USA. 

40.  H.  Wendland,  Scattered  Data  Approximation,  vol.  17.  Cambridge  LIniversity 
Press,  2004.  LInited  Kingdom. 

41.  J.  A.  Reeger  and  B.  Fornberg,  “Numerical  Quadrature  Over  the  Surface  of  a 
Sphere.”  Unpublished  Powerpoint,  2015. 

42.  J.  A.  Reeger,  “Handout:  4  November  2015.”  Private  communication,  Nov.  2015. 

43.  P.  Persson  and  G.  Strang,  “A  Simple  Mesh  Generator  in  MATLAB,”  SIAM  Rev., 
vol.  46,  pp.  329-345,  June  2004. 


82 


44.  J.  A.  Reeger,  “Handout:  29  June  2015.”  Private  communication,  June  2015. 


45.  R.  C.  Yates,  Cassinian  Curves ,  pp.  8-11.  J.  W.  Edwards,  1952.  Ann  Arbor,  MI, 
USA. 

46.  D.  J.  Griffiths  and  R.  College,  The  Three-Dimensional  Delta  Function ,  vol.  3, 
p.  50.  Prentice  Hall  Upper,  1999.  Saddle  River,  NJ,  USA. 

47.  J.  A.  Reeger,  “Handout:  28  January  2016.”  Private  communication,  Jan.  2016. 

48.  B.  Fornberg,  E.  Lehto,  and  C.  Powell,  “Stable  Calculation  of  Gaussian- Based 
RBF-FD  Stencils,”  Comput.  Math.  Appl. ,  vol.  65,  pp.  627-637,  2013. 

49.  B.  Fornberg  and  G.  Wright,  “Stable  Computation  of  Multiquadratic  Interpolants 
for  all  Values  of  the  Shape  Parameter,”  Comput.  Math.  Appl.,  vol.  48,  pp.  853- 
867,  2004. 

50.  B.  Fornberg  and  C.  Piret,  “A  Stable  Algorithm  for  Flat  Radial  Basis  Functions 
on  a  Sphere,”  SIAM  J.  Sci.  Comput.,  vol.  30,  pp.  60-80,  2007. 

51.  R.  Schaback,  Comparisons  of  Radial  Basis  Function  Interpolants ,  ch.  18,  pp.  293- 
305.  World  Scientific,  1993.  Singapore. 

52.  T.  A.  Driscoll  and  B.  Fornberg,  “Interpolation  in  the  Limit  of  Increasingly  Flat 
Radial  Basis  Functions,”  Comput.  Math.  Appl,  vol.  43,  pp.  413-422,  2002. 


83 


REPORT  DOCUMENTATION  PAGE 


Form  Approved 
OMB  No.  0704-0188 


